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Abstract 

In  this  paper  we  develop  a  macroscopic  framework  quantifying  the  hysteresis  and  constitutive 
nonlinearities  inherent  to  ferromagnetic  materials.  In  the  first  step  of  the  development,  we  construct 
Helmholtz  and  Gibbs  energy  relations  at  the  mesoscopic  or  lattice  level  based  on  the  assumption 
that  magnetic  moments  or  spins  are  restricted  to  two  orientations.  Direct  minimization  of  the  Gibbs 
energy  yields  local  average  magnetization  relations  appropriate  for  operating  regimes  in  which  re¬ 
laxation  mechanisms  are  negligible  whereas  the  balance  of  the  Gibbs  and  relative  thermal  energies 
through  Boltzmann  principles  provides  local  models  which  incorporate  mechanisms  such  as  thermal 
after-effects.  To  construct  macroscopic  relations  that  incorporate  material  nonhomogeneities,  poly¬ 
crystallinity,  and  variable  effective  fields,  we  employ  stochastic  homogenization  techniques  based  on 
the  assumption  that  parameters  such  as  local  coercive  and  interaction  fields  are  manifestations  of 
underlying  distributions.  The  resulting  framework  quantifies  in  a  natural  manner  the  anhysteretic 
magnetization  provided  by  decaying  AC  fields  and  guarantees  the  closure  of  biased  minor  loops  once 
transient  accommodation  and  after-effects  are  complete.  Furthermore,  noncongruency  is  achieved 
with  certain  choices  for  the  energy  functionals.  Hence  the  framework  provides  an  energy  basis  for 
certain  extended  Preisach  models  and  the  relation  of  the  framework  to  several  macroscopic  hystere¬ 
sis  models  is  detailed.  The  behavior  of  both  the  nonlinear  anhysteretic  relations  and  full  hysteresis 
model  are  validated  through  comparison  with  experimental  steel  and  nickel  data. 
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1  Introduction 


The  modeling  of  hysteresis  in  ferromagnetic  materials  is  a  classical  problem,  dating  back  to  at  least 
Maxwell  [1],  which  has  profound  implications  for  present  and  projected  applications  utilizing  mag¬ 
netic  materials.  As  the  range  of  magnetic  applications  grows,  so  too  does  the  list  of  requirements 
necessary  for  the  models  used  to  quantify  hysteresis  and  constitutive  nonlinearities  inherent  to  the 
compounds.  For  example,  transducer  design,  material  characterization,  optimization  of  magnetic 
recording  media,  and  development  of  magnetic  computing  paradigms  require  high  fidelity  models 
that  are  feasible  to  implement  in  optimization  algorithms.  Real-time  implementation  of  model-based 
control  algorithms  for  such  systems  necessitates  the  development  of  low-order  macroscopic  models 
which  are  easily  constructed  and  updated  to  accommodate  changing  environmental  conditions.  To 
optimize  material  design,  however,  it  is  necessary  to  construct  models  at  the  microscopic  and  meso¬ 
scopic  levels  to  quantify  and  predict  the  effects  of  changing  material  composition.  Moreover,  these 
microscopic  models  must  be  commensurate  with  macroscopic  measurements  to  permit  evaluation 
and  updating  of  the  material  designs.  This  necessitates  the  development  of  multiscale  modeling 
hierarchies  in  which  energy-based  fine-scale  models  are  used  to  predict  effective  parameters  in  higher 
level  models. 

As  a  prelude  to  constructing  any  energy-based  models,  it  is  necessary  to  at  least  qualitatively 
understand  the  physical  mechanisms  which  produce  hysteresis  and  constitutive  nonlinearities.  In 
magnetic  materials,  these  properties  are  due  to  a  number  of  mechanisms  but  the  primary  sources  of 
hysteresis  are  moment  or  domain  interactions,  domain  wall  losses  caused  by  material  inclusions,  and 
material  anisotropies  [2-5] .  Micromagnetic  models  can  accommodate  a  number  of  these  mechanisms 
but  do  so  at  high  computational  cost.  There  are  presently  no  macroscopic  models  which  quantify 
all  of  these  effects  for  general,  broadband,  variable  temperature,  variable  stress  operating  conditions, 
and  present  macroscopic  theories  address  one  or  two  of  these  mechanisms  as  dictated  by  the  regimes 
under  consideration. 

In  this  paper,  we  develop  a  macroscopic  framework  which  quantifies  hysteresis  and  constitutive 
nonlinearities  in  the  H-M  and  H-B  relation  for  a  variety  of  ferromagnetic  materials  and  yields  an- 
hysteretic  magnetization  relations  in  a  natural  manner.  It  guarantees  the  closure  of  biased  minor 
loops  once  after-effect  and  accommodation  processes  are  complete  but  does  not  enforce  congruency  - 
in  accordance  with  noncongruencies  measured  in  certain  operating  regimes.  The  model  incorporates 
relaxation  mechanisms  but  neglects  eddy  current  losses  so  it  should  be  employed  for  low  frequency 
material  characterization  or  architectures  for  which  eddy  current  losses  are  minimal  (e.g.,  laminated 
Terfenol-D  rods).  It  is  constructed  in  the  context  of  uniaxial  moment  configurations  but  does  not 
incorporate  additional  anisotropy  energy  mechanisms;  hence  it  applies  to  isotropic  polycrystalline 
materials,  a  variety  of  uniaxial  regimes,  and  weakly  anisotropic  materials.  Finally,  the  characteriza¬ 
tion  framework  is  sufficiently  efficient  to  permit  model-based  system  and  control  design. 

In  the  first  step  of  the  model  development,  a  mean  field  approximation  for  the  exchange  energy 
of  moments  on  a  uniform  lattice  is  balanced  with  entropy  effects  to  construct  a  Helmholtz  free 
energy  relation.  This  in  turn  is  used  to  construct  a  piecewise  quadratic  Helmholtz  energy  which 
retains  salient  features  of  the  statistical  mechanics  model  at  fixed  temperatures  but  is  more  efficient 
to  implement.  The  incorporation  of  the  magnetoelastic  energy  subsequently  yields  a  Gibbs  energy 
relation  which  quantifies  changes  in  the  energy  due  to  an  applied  field.  For  general  operating  regimes, 
the  Gibbs  and  relative  thermal  energies  are  balanced  through  Boltzmann  probability  relations  to 
yield  a  local  average  magnetization  model  which  incorporates  the  thermal  activation  mechanisms 
which  produce  thermal  after-effects  [2]  and  certain  manifestations  of  accommodation  or  reptation  [6]. 
In  the  limit  of  negligible  thermal  activation,  it  is  proven  that  these  relations  reduce  to  the  local 
magnetization  kernels,  or  hysterons,  obtained  by  minimizing  the  Gibbs  energy.  For  certain  regimes, 
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it  is  illustrated  that  the  kernels  are  given  by  Ising  relations.  In  the  second  step  of  the  development, 
we  incorporate  lattice  nonhomogeneities,  inclusions,  polycrystallinity,  texture  and  certain  mean  field 
effects  by  assuming  that  quantities  such  as  the  local  coercive  and  effective  fields  are  manifestations 
of  underlying  distributions  rather  than  fixed  parameters.  Stochastic  homogenization  in  this  manner 
yields  macroscopic  models  suitable  for  material  characterization. 

To  place  the  modeling  framework  in  perspective,  we  compare  it  with  four  presently  employed 
macroscopic  models:  Jiles-Atherton,  Stoner- Wohlfarth,  Globus  and  Preisach  —  detailed  comparison 
of  these  models  can  be  found  in  [7].  These  four  approaches  were  chosen  to  illustrate  similarities 
and  differences  with  the  proposed  framework,  and  they  represent  what  can  now  be  considered  as 
established  macroscopic  hysteresis  models.  However,  they  are  by  no  means  exhaustive  and  additional 
hysteresis  models  can  be  found  in  [2-4, 8,9].  Additionally,  details  about  micromagnetic  theory  and 
its  relation  to  the  proposed  model  is  included  in  Section  2  where  it  is  more  in  context. 

Jiles-Atherton  Model 

The  Jiles-Atherton  model  characterizes  hysteresis  in  isotropic  materials  through  the  quantifica¬ 
tion  of  domain  wall  losses  [4, 10-12].  In  its  original  formulation,  the  model  was  constructed  in  two 
steps:  (i)  quantification  of  the  anhysteretic  magnetization,  and  (ii)  incorporation  of  irreversible  and 
reversible  domain  wall  effects.  The  anhysteretic  magnetization  is  modeled  through  the  Langevin 
relation 

Man  =  Ms  [coth {He/a)  -  ( a/He)\ 

He  =  H  +  aM  ^ 

where  Ms  is  the  saturation  magnetization,  cc  is  a  mean  field  parameter,  and  a  is  a  parameter  having 
dimensions  of  field.  The  irreversible  and  reversible  domain  wall  losses  are  incorporated  by  determining 
the  magnetostatic  energy  required  to  reorient  moments  in  the  presence  of  the  effective  field  He.  In 
this  manner,  the  model  also  incorporates  certain  rotational  effects. 

The  original  model  has  subsequently  been  extended  to  include  eddy  current  losses  [13],  certain 
anisotropies  [14]  and  to  enforce  closure  of  biased  minor  loops.  The  minor  loop  extension  provided 
by  Carpenter  [15]  is  phenomenological  in  the  sense  that  it  is  based  on  translates  and  scaling  of  the 
major  loop  curves  whereas  the  extension  provided  by  Jiles  [16]  relies  on  a  priori  knowledge  of  future 
turning  points.  This  reduces  the  utility  of  the  model  for  dynamic  control  applications  where  turning 
points  are  determined  by  a  feedback  law  responding  to  state  measurements  and  hence  are  unknown 
before  the  control  commences. 

It  will  be  observed  that  the  proposed  model  yields  a  hysteresis  kernel  at  the  mesoscopic  level 
which  is  formulated  as  the  Ising  relation  M  =  Mst&nh(He/a)  and  hence  agrees  through  first-order 
terms  with  (1).  The  primary  difference  between  the  theories  lies  in  underlying  energy  formulation 
and  manner  through  which  losses  are  incorporated. 

Stoner-Wohlfarth  Model 

The  original  Stoner-Wohlfarth  model  quantified  the  rotation  of  noninteracting,  single-domain 
particles  having  uniaxial  anisotropy  [17].  While  the  model  has  received  widespread  use  in  the  mag¬ 
netic  recording  industry,  its  use  for  general  material  characterization  was  originally  limited  by  the 
fact  that  it  did  not  incorporate  moment  interactions.  A  number  of  recent  extensions  to  both  the 
model  and  underlying  philosophy  have  substantially  improved  its  utility.  The  model  was  extended 
to  include  cubic  anisotropies  by  Lee  and  Bishop  [18]  whereas  Armstrong  [19],  Clark  et  al.  [20],  and 
Jiles  and  Thoelke  [21]  extended  the  model  to  quantify  magnetoelastic  effects  in  Terfenol-D.  Certain 
mean  field  effects  are  incorporated  in  [22,23]  while  pinning  losses  are  incorporated  in  [7]  where  it  is 
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illustrated  that  this  latter  mechanism  is  necessary  to  achieve  physical  minor  loop  behavior.  Finally, 
relations  between  the  Stoner- Wohlfarth  model  and  micromagnetic  models  are  detailed  in  [8]. 

The  hysteresis  kernels,  or  hysterons,  provided  by  the  proposed  model  are  analogous  to  those 
predicted  by  the  Stoner-Wohlfarth  model  when  the  applied  magnetic  field  is  aligned  with  the  easy 
axis  of  particles.  This  is  consistent  with  the  assumption  in  the  proposed  framework  that  moment 
alignment  occurs  in  two  diametrically  opposite  directions.  The  assumptions  differ  from  the  original 
Stoner-Wohlfarth  model  in  that  moment  interactions  are  incorporated  whereas  anisotropy  energies 
are  neglected. 

Globus  Model 

The  Globus  model  preceded  that  of  Jiles  of  Atherton  and  bears  certain  similarities  in  that  it 
quantifies  irreversible  and  reversible  hysteresis  effects  through  domain  wall  mechanisms  [24],  A 
primary  difference  between  the  theories  lies  in  the  Globus  assumption  that  domain  walls  are  pinned 
on  grain  boundaries  which  results  in  the  simplifying  assumption  that  reversible  domain  wall  bending 
and  irreversible  domain  wall  displacement  can  be  considered  on  a  single,  representative  spherical 
grain.  While  efficient  to  implement,  this  model  lacks  a  number  of  mechanisms  required  for  high 
performance  applications. 

The  anhysteretic  curves  provided  by  the  Globus  model  are  analogous  to  the  mesoscopic  hysteresis 
kernel  provided  by  the  new  theory.  However,  the  latter  includes  moment  interactions  and  physical 
minor  loop  behavior  which  makes  it  advantageous  for  general  applications. 


Preisach  Framework 

Preisach’s  original  model  quantified  the  hysteretic  relation  between  input  fields  H{t)  and  the 
magnetization  M(t)  through  a  superposition  of  rectangular  hysteresis  relays  or  kernels  [ksljS2(H)](t) 
[25].  Here  si  and  S2,  with  .si  <  s 2,  provide  thresholds  at  which  the  kernel  switches  between  +1  and 
— 1.  The  magnetization  is  modeled  as 

POO  POO 

[M(H)](t)=  /  /  u(r,s)[ks-rtS+r(H)](t)dsdr  (2) 

J  0  J  —  00 


where  v  is  a  density,  or  weighting  function  which  depends  on  properties  of  the  material  under 
consideration  [26].  It  is  illustrated  in  [7,9]  that  one  choice  of  v  is  the  normal  density  function 


v(h!,hc)  =  Ms  e-{hc-hc)2/2a2ce-hj/2aj 


27T  (J  q(J  j 


(3) 


where  hi  and  hc  denote  interaction  and  coercive  fields,  hc  is  an  average  coercive  field,  and  07,  ac  are 
respective  variances.  For  general  characterization,  the  weight  or  measure  v  can  be  estimated  from 
measured  data  through  techniques  analogous  to  those  described  in  [27,28]. 

The  advantage  of  the  Preisach  methodology  lies  in  its  generality.  From  a  solely  mathematical 
perspective,  it  can  be  used  to  characterize  material  behavior  in  regimes  where  the  underlying  physics 
is  poorly  understood  or  unknown.  It  also  yields  models  which  can,  at  least  approximately,  be  in¬ 
verted  for  linear  control  design  [28].  However,  these  advantages  are  complemented  by  a  number  of 
disadvantages.  First,  the  nonphysical  nature  of  the  measures  or  parameters  makes  it  difficult  to 
construct  models  using  known  physical  behavior  or  to  employ  attributes  of  the  data  for  updating 
models  to  accommodate  changing  operating  conditions.  While  extensions  to  the  Preisach  model 
have  recently  been  proposed  to  facilitate  parameter  identification  through  correlation  with  physical 
mechanisms  [9,  29] ,  the  accurate  characterization  of  biased  minor  loops  typically  requires  general 
weights  comprised  of  a  large  number  of  nonphysical  parameters.  A  more  fundamental  difficulty 
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arises  in  operating  regimes  involving  broadband  input  signals  and  temperature  or  load  dependence 
since  the  classical  Preisach  model  is  based  on  the  assumptions  of  frequency,  load  and  temperature 
invariance.  Finally,  classical  Preisach  models  erroneously  enforce  congruency  in  certain  regimes  and 
do  not  accommodate  reversible  magnetization  mechanisms.  As  detailed  in  [6,7,9,30],  extensions  to 
the  classical  theory  have  been  developed  to  address  a  number  of  these  issues.  However,  these  exten¬ 
sions  add  significant  complexity  to  the  theory  and  reduce  its  efficiency  for  real-time  implementation. 
For  example,  one  technique  for  incorporating  temperature-dependence  is  to  employ  vector-valued 
measures  or  weights  [v(r,  s)](T)  (e.g.,  see  [31]).  For  control  design,  however,  this  necessitates  the  use 
of  lookup  tables  which  reduces  significantly  the  efficiency  of  associated  control  algorithms. 

As  illustrated  in  [32,33],  the  proposed  model  can  be  interpreted  as  providing  an  energy  basis  for 
certain  extended  Preisach  models.  In  summary,  the  kernels  derived  through  energy  considerations 
at  the  mesoscopic  level  provide  the  Preisach  hysterons  whereas  the  stochastic  densities,  used  to 
obtain  macroscopic  models,  provide  the  Preisach  weights.  However,  the  proposed  framework  differs 
from  the  classical  Preisach  model  in  five  crucial  aspects,  (i)  The  first  is  the  energy  basis  of  the 
model  —  formulation  through  energy  analysis  provides  a  low-order  model  which,  for  certain  density 
choices,  has  parameters  that  can  be  physically  correlated  with  properties  of  the  data,  (ii)  Due  to 
the  energy  basis  of  the  framework,  stress  and  temperature-dependencies  are  incorporated  in  the 
basis  or  kernel  in  the  new  model  (e.g.,  see  [34,35])  whereas  they  enter  the  weights  in  the  Preisach 
formulation.  Because  they  are  incorporated  in  the  kernel,  the  model  automatically  incorporates 
these  effects  which  eliminates  the  necessity  of  vector-valued  weights  or  lookup  tables.  From  the 
perspective  of  implementation,  this  indicates  that  only  one  set  of  parameters  must  be  identified 
for  the  proposed  model  and  no  switching  between  parameter  sets  is  required  during  operation.  This 
significantly  augments  the  efficiency  of  characterization  and  control  algorithms  employing  the  model, 
(iii)  The  incorporation  of  relaxation  mechanisms  through  the  energy  basis  provides  the  framework 
with  the  property  that  it  accommodates  nonclosure  or  nondeletion  properties  in  accordance  with 
measured  material  properties.  In  this  case,  the  framework  provides  an  energy  basis  for  certain 
extended  Preisach  formulations  based  on  Arrhenius  relations  [9].  (iv)  Derivation  of  the  theory  from 
Boltzmann  principles  yields  kernels  or  hysterons  which  accommodate  the  noncongruency  observed 
for  certain  materials  and  operating  regimes.  This  is  in  contrast  to  input  and  output-dependent 
(moving)  Preisach  formulations  which  incorporate  noncongruency  through  input  or  output-dependent 
densities,  (v)  The  model  automatically  incorporates  reversible  magnetization  mechanisms  for  small 
AC  field  excursions  about  a  fixed  DC  value  -  hence  it  accurately  characterizes  reversible  material 
behavior  following  field  reversal  without  the  extensions  required  for  Preisach  formulations. 

The  proposed  framework  builds  upon  and  significantly  extends  ferromagnetic  theory  originally 
developed  in  [36].  The  extensions  include  the  derivation  of  energy  relations  based  on  statistical 
mechanics  tenets,  formulation  of  the  model  in  terms  of  general  densities,  and  rigorous  analysis 
establishing  the  convergence  of  models  incorporating  thermal  activation  to  those  derived  through 
minimization  of  the  Gibbs  energy  as  relative  thermal  energies  become  negligible.  In  the  present 
work,  we  also  establish  the  manner  through  which  the  anhysteretic  magnetization  is  characterized 
and  summarize  highly  efficient  implementation  algorithms. 

Appropriate  Helmholtz  and  Gibbs  energy  relations  are  constructed  at  the  lattice  level  in  Sec¬ 
tion  2  and  used  in  Section  3  to  quantify  the  local  average  magnetization  M  and  local  anhysteretic 
magnetization  Man  in  the  presence  and  absence  of  thermal  activation.  These  local  relations  are 
combined  with  stochastic  homogenization  techniques  in  Section  4  to  construct  macroscopic  models 
which  quantify  the  hysteresis  and  constitutive  nonlinearities  inherent  to  a  variety  of  ferromagnetic 
compounds.  The  properties  of  both  the  anhysteretic  and  full  hysteresis  models  are  illustrated  in 
Section  5  through  numerical  simulations  and  comparison  with  experimental  steel  data. 
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2  Energy  Relations  for  Homogeneous  Materials 


The  microscopic  and  macroscopic  behavior  of  ferromagnetic  materials  is  defined  by  the  exchange, 
anisotropy,  magnetostatic  and  magnetoelastic  energies,  and  hence  a  complete  energy-based  theory 
quantifying  hysteresis  and  magnetic  properties  of  these  materials  must  include,  or  at  least,  accom¬ 
modate  these  contributions.  The  magnetostatic  energy  quantifies  long  range  interactions  between 
magnetic  moments  and  applied  fields  and  is  significant  in  both  ferromagnetic  and  paramagnetic 
materials.  Additionally,  ferromagnetic  materials  exhibit  exchange  interactions  between  neighbor¬ 
ing  atomic  spins  which  serves  to  align  moments.  This  effect  is  short  range  and  typically  influences 
nearest,  or  next-nearest,  neighbors.  The  exchange  interactions  and  associated  exchange  energy  also 
differentiate  ferromagnetic  and  paramagnetic  compounds.  The  anisotropy  energy  quantifies  changes 
in  the  internal  energy  of  materials  due  to  changes  in  the  direction  of  the  magnetization.  Whereas 
magnetic  anisotropy  can  be  due  to  a  number  of  factors,  we  will  focus  primarily  on  crystalline  and 
stress  anisotropies.  Finally,  the  magnetoelastic  energy  quantifies  magnetomechanical  coupling  inher¬ 
ent  to  the  materials. 

The  first  statistical  mechanics  model  for  the  exchange  interactions  was  posed  by  Ising  in  1925 
and  was  based  on  the  assumption  of  a  linear  lattice  of  magnetic  moments  in  which  only  neighbors 
interact  [37].  This  was  later  extended  by  Heisenberg  in  1928  to  include  quantum  effects  and  complete 
correlation  between  neighboring  electrons  having  overlapping  wave  functions  [38] .  In  the  Heisenberg 
model,  the  interaction  energy  between  spins  S,;  and  S;  is 

U  =  -2ySj:  •  S j  (4) 

where  J  is  an  exchange  integral  with  J  >  0  for  ferromagnetic  materials  and  J  <  0  for  antiferro¬ 
magnetic  compounds.  The  exchange  energy  for  the  system  is  then  obtained  by  summing  over  all 
magnetic  moments  which  yields 

U=  -2j2JijSi  -Sj,  (5) 

ij 

The  Ising  model  can  be  obtained  from  that  of  Heisenberg  by  truncating  the  interaction  energy  for 
the  lattice.  For  example,  if  quantization  is  assumed  to  take  place  only  in  the  z-direction  (e.g.,  the 
direction  of  an  applied  field  H),  the  full  inner  product  Sj  •  S j  =  SixSjx  +  SiySjy  +  SizSjz  is  replaced 
in  the  Ising  model  by  the  z-component  SizSjz.  If  the  restricted  spins  are  denoted  by  cr*  =  ±1,  the 
Ising  relation  for  the  exchange  energy  can  be  formulated  as 

U  =  —2  ^2  Jij&iVj  (6) 

(ij) 

where  the  notation  (ij)  indicates  summation  over  nearest  neighbors.  In  addition  to  the  simplification 
provided  by  reduced  dimensionality,  the  Ising  model  admits  a  classical  treatment  of  the  spins,  due  in 
part  to  the  fact  that  spins  commute  in  the  truncated  expansion  and  hence  can  be  treated  as  variables, 
whereas  spins  must  be  treated  as  quantum  mechanical  operators  in  the  Heisenberg  model.  While 
the  Ising  model  proves  sufficiently  accurate  when  characterizing  the  exchange  energy  in  a  number  of 
applications,  it  is  necessary  to  include  the  quantum  effects  incorporated  in  the  Heisenberg  model  to 
quantify  mean  field  properties  from  first  principles. 

As  detailed  in  [8,39],  the  Heisenberg  energy  relation  (5)  is  exactly  solvable  for  only  a  few  cases, 
one  of  which  is  the  Ising  model,  and  is  completely  isotropic.  The  inclusion  of  the  anisotropy  and 
magnetostatic  energy  in  the  Heisenberg  Hamiltonian  or  energy  formulation  yields  a  model  that  is 
prohibitively  expensive  to  approximate  due  to  its  generality.  This  necessitates  the  consideration  of 
simplified  theories  motivated  by  the  quantum  energy  relations  but  tractable  for  implementation. 
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One  technique  in  this  vein  is  the  theory  of  micromagnetics  which  originated  with  the  work  by 
Landau  and  Lifshitz  in  1935  [40]  on  the  analysis  of  domain  walls.  Significant  contributions  to  the 
theory  were  provided  by  Brown  [41-43]  and  subsequent  researchers;  e.g.  see  [8].  The  basic  tenet  in 
this  approach  is  to  employ  classical  theory  in  combination  with  quantum  principles  to  predict  the 
distribution  of  spins  through  the  minimization  of  general  energy  relations 

U  =  UE  +  UA  +  UX  +  Um  (7) 

where  Ue,  Ua,  U\  and  Um  respectively  denote  the  exchange,  anisotropy,  magnetoelastic  and  magneto¬ 
static  energies.  The  capability  this  theory  provides  for  predicting  domain  structures  in  ferromagnetic 
materials  is  illustrated  in  [44]  where  a  micromagnetic  energy  relation  of  the  form  (7),  with  the  specific 
energy  components  defined  by 

x  y  z 

UA  =  —2  EXE  [ill  (afa%  +  +  a|a^)  +  K20l\ol\o^ 

x  y  z 

o  '  (8) 

Ux  =  2  </> 

x  y  z 

x  y  z 

was  minimized  on  a  Cray  X-MP/22  for  a  cubic  grid  consisting  of  22  x  22  x  22  exchange  coupled 
spins.  Here  K\  and  are  anisotropy  constants  and  an,  a.2  and  a 3  are  direction  cosines  of  a  moment 
m.j  at  the  center  of  the  cube.  Furthermore,  A,  cr  and  </>  respectively  denote  the  magnetostriction 
constant,  magnetoelastic  stress  and  angle  between  the  magnetization  and  stress.  Finally,  Hj  denotes 
the  interaction  field  at  m;  due  to  all  moment. 

In  addition  to  the  capabilities  that  this  theory  provides  for  predicting  the  domain  structure  and 
dynamics  of  ferromagnetic  materials,  it  is  established  in  [8]  that  this  approach  provides  a  set  of 
differential  equations  for  which  the  Stoner- Wohlfarth  model  is  an  eigenfunction  or  mode. 

However,  due  to  their  complexity,  models  derived  through  micromagnetic  principles  presently 
preclude  real-time  implementation.  The  predictive  capabilities  provided  by  the  models  include  more 
detail  than  is  necessary  for  modeling  hysteresis  in  a  manner  that  facilitates  system  and  control  design. 
To  provide  such  quasi- macroscopic  models,  we  construct  an  energy  formulation  which  incorporates 
certain  aspects  of  the  Ising  and  micromagnetics  models  at  the  microscopic  scale  but  employs  statisti¬ 
cal  techniques  rather  than  physical  principles  to  obtain  a  commensurate  macroscopic  magnetization 
model. 

2.1  Helmholtz  Energy 

We  consider  two  techniques  for  specifying  the  Helmholtz  energy  at  the  lattice  level.  The  first  is  based 
on  statistical  mechanics  principles  and  hence  includes  certain  temperature  dependencies;  the  second 
is  obtained  through  an  approximation  of  the  first  for  fixed  temperature  regimes. 

The  statistical  mechanics  model  is  based  on  an  approximation  to  the  Ising  model  which  has  the 
requisite  assumption  that  spins  or  magnetic  moments  are  restricted  to  two  possible  orientations, 

(Ti  =  ±1. 

The  assumption  that  spins  have  two  preferred  orientations  appears  at  first  to  be  highly  restrictive 
but  can  in  fact  be  physically  motivated  for  a  number  of  regimes.  From  a  classical  perspective, 
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Figure  1:  (a)  Crystallographic  orientation  of  magnetic  moments  in  cobalt  and  iron  (after  [4]).  (b)  Cor¬ 
responding  domain  structure  which  reflects  the  uniaxial  or  hexagonal  anisotropy  of  cobalt  and  the 
cubic  anisotropy  of  iron. 


this  will  be  at  least  approximately  true  for  materials  having  uniaxial  crystalline  anisotropies  or 
systems  in  which  uniaxial  stresses  dominate  the  crystalline  structure.  Materials  exhibiting  uniaxial 
crystalline  anisotropies  include  cobalt  and  a  number  of  rare  earth  metals  and  alloys  (e.g.,  Terbium 
single  crystals).  This  produces  domain  structures  in  which  moments  are  highly  parallel  or  antiparallel 
as  depicted  in  Figure  1. 

To  illustrate  a  regime  in  which  stresses  dominate  crystalline  anisotropies,  consider  the  Terfenol- 
D  transducer  depicted  in  Figure  2(a)  and  detailed  in  [45-47].  In  present  manufacturing  processes, 
Terfenol-D  crystals  are  grown  in  Dendrite  sheets  oriented  in  the  [112]  directions  as  depicted  in 
Figure  2(b).  At  the  prestress  levels  employed  in  present  transducer  designs,  the  preferred  orientation 
of  domains  is  shifted  from  the  original  eight  (111)  magnetic  easy  axes  to  the  two  axes  [111]  and  [111] 
perpendicular  to  the  [112]  axis  of  the  rod.  In  the  presence  of  a  field  H  generated  by  an  applied  current 
/  to  the  solenoid,  moments  first  rotate  irreversibly  to  the  [111]  easy  axis  and  then  rotate  reversibly 
to  the  [112]  axis.  For  these  transducer  constructions,  stress  anisotropies  can  dominate  crystalline 
anisotropies  to  provide  regimes  for  which  the  assumption  of  two  spin  orientations  provides  reasonable 
approximations. 


Wound  Wire  Solenoid 


Compression 
Bolt - - 


Permanent  Magnet 


Figure  2:  (a)  Cross  section  of  a  typical  Terfenol-D  magnetostrictive  transducer,  (b)  Orientation  of 
Terfenol-D  crystals. 


7 


Finally,  this  assumption  can  be  motivated  by  noting  that  from  a  quantum  perspective,  spins 
cannot  be  uniformly  oriented  and  one  allowed  orientation  is  parallel  and  opposite  to  an  applied  field. 
This  interpretation  can  be  used  to  explain  the  accuracy  of  the  theory  for  quantifying  the  behavior  of 
certain  materials  such  as  iron  which  has  the  crystalline  and  domain  properties  illustrated  in  Figure  1 
and  exhibits  cubic  anisotropy. 


2.1.1  Temperature-Dependent  Helmholtz  Energy 

The  statistical  mechanics  model  is  based  on  an  approximation  to  the  Ising  model  first  proposed 
by  Gorsky  in  the  analysis  of  order-disorder  transitions  in  binary  alloys  [48].  In  this  context,  it 
was  later  extended  by  Bragg  and  Williams  to  include  the  concept  of  long  range  order  [49-51].  An 
underlying  tenet  in  the  Bragg-Williams  theory,  which  simplifies  subsequent  computations,  is  the 
assumption  that  the  energy  of  an  individual  atom  is  determined  by  the  average  order  of  the  system 
rather  than  the  fluctuating  states  of  adjacent  atoms.  For  this  reason,  the  model  is  often  termed 
the  mean  field  or  molecular  field  approximation  to  the  Ising  model.  To  construct  a  ferromagnetic 
model,  we  make  the  same  assumption  regarding  magnetic  moments  or  spins.  Further  details  about 
this  approach,  including  some  discussion  concerning  its  application  to  ferromagnetic  materials,  can 
be  found  in  [52, 53]. 

We  consider  an  arbitrary  lattice  of  volume  V  and  mass  v  comprised  of  N  =  N+  +  N_  cells,  each  of 
which  is  assumed  to  contain  one  spin  or  magnetic  moment.  In  accordance  with  the  Ising  assumptions, 
the  spin  orientations  are  constrained  to  be  cq  =  ±1,  and  N+  and  7V_  respectively  denote  the  number 
of  positive  and  negative  spins  in  the  lattice.  We  note  that  due  to  the  initial  assumption  of  material 
homogeneity,  this  lattice  structure  is  representative  of  that  found  throughout  the  structure.  If  each 
spin  has  a  moment  m,  the  magnetization  for  the  lattice  is 


M 


m 

V 


N 


i=  1 


N+-N. ), 


from  which  it  follows  that 


N+  = 


N 

2 


(9) 


(10) 


Here  Ms  =  Nm/V  denotes  the  technical  saturation  magnetization  which  occurs  when  all  moments 
are  aligned.  Additionally,  we  make  the  assumption  that  only  adjacent  moments  interact. 

To  quantify  the  energy  required  to  reorient  moments,  we  employ  the  mean  field  approximation 
of  Bragg  and  Williams  and  make  the  assumption  that  the  average  exchange  energy  is  proportional 
to  M/Ms ;  that  is, 

$  =  $0  M/Ms  (11) 

where  $o  denotes  the  energy  required  to  reorient  a  single  moment  if  the  lattice  is  completely  ordered 
(M  =  Ms).  For  the  case  of  a  homogeneous  lattice,  $0  is  considered  to  be  constant.  For  nonho- 
mogeneous  and  polycrystalline  materials,  $0  will  be  considered  as  a  manifestation  of  an  underlying 
statistical  distribution  as  discussed  in  Section  4.  We  also  note  that  $0  is  related  to  the  exchange 
integral  J  employed  in  (6)  through  the  expression 


$0  =  2  tJ 


(12) 


where  £  denotes  the  number  of  neighbors  adjacent  to  a  site.  Hence  £  =  2  for  a  1-D  lattice  chain, 
£  =  4  for  a  2-D  rectangular  lattice,  and  £  =  6,8  or  12,  respectively,  for  3-D  cubic,  body-centered 
cubic,  or  face-centered  cubic  lattices.  The  fact  that  (11)  is  independent  of  the  exact  lattice  structure 
has  led  Pathria  to  refer  to  the  subsequent  model  as  a  zeroth  approximation  of  the  Ising  model  [53] . 

We  now  consider  the  decrease  in  internal  energy  due  to  a  change  from  N+  to  N+  +  dN+.  From 
the  mean  field  approximation  (11),  each  switch  requires  ^ojj~  in  energy  so  the  change  in  internal 
energy  for  a  unit  volume  is 


dUE 


$0  M 
VM~S 


dN . 


<t>oM 

VMS 


N 

2  Ms 


dM 


(13) 


where  the  second  equality  follows  from  (10).  Integration,  in  combination  with  (10),  yields  the  relation 


UE 


$o N  (  M2\ 
4V  V  Mi) 


+  U0 


(14) 


for  the  exchange  energy.  Since  we  are  interested  in  relative  rather  than  absolute  measures  of  energy, 
we  take  Uq  =  0  which  specifies  that  the  completely  ordered  state  has  an  internal  energy  of  zero. 

As  detailed  in  [52],  the  entropy  S  for  the  system  is  given  by 

kN 

S=yr  kl  W  (15) 

where  k  is  Boltzmann’s  constant  and  W  quantifies  the  number  of  ways  moments  can  be  arranged  in 
the  lattice  to  yield  the  magnetization  M.  By  noting  that  this  is  equivalent  to  arranging  N+  moments 
in  N  sites,  and  employing  Stirling’s  approximation 


In  x!  =  x  In  x  —  x  . 


(16) 


the  entropy  can  be  formulated  as 


5  - 


In 


V 

kN 

V 


N 

N+ 

N\ 


N-\  N+\ 


ln2-l  +  M/M  /  M 

2  V  M. 


1  —  M/Ms  /  M 
2  V  M. 


+  S0 


(17) 


-kN 
2 VAN 


,  ,M  +  MS\  [  f 

M  In  (  - w  1  +  Ms  In  (  1  -  f  —  J 


Ms  -  M 


+  So 


where  So  =  In  2. 

The  Helmholz  energy  for  the  lattice  is  then 
ip(M,  T)  =  U  -  ST 


[1  -  (M/Ms)2]  + 


$0 N  . . on  TkN 

4V 

HhM. 


[1  -  (M/Ms)2]  + 


2  VMS 
HhT 


Min  (  M  +  Mt  ]  +  Ms  In  (l  -  (M/Ms)2) 


2  Tr 


M  In 


Ms  -  M 

M  +  Ms 
Ms  -  M 


+  Ms  In  (l  —  (M /Ms)2) 


(18) 
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(a)  (b) 

Figure  3:  Helmholtz  energy  specified  by  (18)  for  (a)  T  <  Tc,  and  (b)  T  >  Tc. 


where  H ^  =  .y^y  is  a  bias  field  and  Tc  =  ^  denotes  the  Curie  temperature.  The  initial  assumption 
that  the  exchange  energy  $o  is  constant  implies  that  Hj,  will  also  be  constant  for  homogeneous 
materials.  This  assumption  will  be  relaxed  in  Section  4  to  include  statistically  distributed  values  of 
Hh  for  modeling  nonhomogeneous  and  polycrystalline  materials. 

As  illustrated  in  Figure  3,  the  Helmholtz  relation  (18)  exhibits  double  well  behavior  for  temper¬ 
atures  T  <TC  and  single  well  behavior  for  T  >TC.  This  is  consistent  with  the  transition  exhibited 
between  ferromagnetic  and  paramagnetic  phases. 


2.1.2  Temperature-Invariant  Helmholtz  Energy 

Whereas  the  Helmholtz  relation  (18)  incorporates  a  number  of  the  properties  desired  for  microscopic 
material  characterization,  the  logarithmic  components  add  complexity  to  resulting  macroscopic  mod¬ 
els  which  can  reduce  the  efficiency  of  algorithms  when  considered  for  real-time  implementation.  For 
applications  requiring  high  efficiency,  a  simplified  Helmholtz  relation  can  be  obtained  by  retaining 
the  quadratic  behavior  of  (18)  for  fixed  temperature  regimes. 

To  determine  an  appropriate  piecewise  quadratic  model,  we  consider  the  Taylor  expansion 

4 >(M,T)  =  ip(Mo,  T)  +  (M  -  M°)|^(Mo,T)  +  (M  ~M°)2  (Mq>  t)  +  0(M3)  (19) 

for  fixed  T  <  Tc,  where  Mq  is  taken  to  be  an  equilibrium,  and 


dip 

dM 


~ Hh 
Ms 


M  + 


~wr~tanh  l(M/Ms) , 

J-  C. 


d2ip  _  -Hh  HhT 

dM 2  Ms  +  TCMS  [1  -  (M/Ms)2]  ' 


(20) 


From  the  necessary  condition 
solutions  to 


,T) 


0,  the  equilibria  are  determined  to  be  the  two  stable 


M  =  Ms  tanh 


(21) 
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in  addition  to  the  unstable  solution  M  =  0.  The  parameters  a  and  a(T )  are  specified  by 


Hh  frn  HhT 

M  ’  K  ’  Tr 


(22) 


If  we  let  Mr  and  —Mr  denote  locations  of  the  stable  equilibria  determined  through  solution  of  (21), 
as  depicted  in  Figure  3(a),  then  it  can  be  directly  established  that  the  quadratic  approximations  to 
(19)  in  neighborhoods  of  the  equilibria  Mq  =  0,  —Mr  and  Mr  are 


where 


ip(M,  T)  = 


Cl 


'  ci  —  kf(T)M2 
<  c2(T)  +  k22(T)(M  +  MR)2 

k  c2(T)  +  k2(T)(M  -  Mr)2 

=  HhMs  b2(T)  =  Hh 

2  ’  '  2 MSTC 


,  Mo  =  0 


,  M0  =  -Mr 


,  M0  =  Mr 


(Tc  -  T) 


(23) 


(24) 


with  analogous,  but  more  complicated,  expressions  for  c2{T)  and  k2(T). 

For  fixed  temperature  regimes,  this  motivates  the  consideration  of  the  piecewise  quadratic  defi¬ 
nition 


-ip(M) 


r  \r](M  +  Mr)2 
J  \r]{M  —  Mr)2 


,  M  <  -Mi 
,  M  >  M/ 


(25) 


[  -  Mr)  -  Mr)  ,  |M|  <  Ml 


as  a  second  choice  for  the  Helmholtz  energy.  As  illustrated  in  Figure  4(a),  Mj  and  Mr  respectively 
denote  the  inflection  point  and  magnetization  at  which  the  minimum  of  ip  occurs.  It  will  be  estab¬ 
lished  in  subsequent  discussion  that  Mj  and  Mr  represent  parameters  to  be  estimated  through  a 
least  squares  fit  to  data  when  quantifying  specific  materials. 


V(M)=G(0,M) 
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Figure  4:  (a)  Helmholtz  energy  ?/>  and  Gibbs  energy  G  for  increasing  field  H  ( H2  >  H\  >  0). 
(b)  Dependence  of  the  local  average  magnetization  M  on  the  field  in  the  absence  of  thermal  activation. 
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2.2  Gibbs  Energy 

The  Helmholtz  relations  (18)  or  (25)  quantify  certain  aspects  of  the  exchange  energy  Ue  for  ferro¬ 
magnetic  materials.  To  incorporate  the  work  done  by  an  applied  field,  we  note  from  (8)  that  the 
magnetostatic  energy  can  be  expressed  as  Um  =  hom  •  H,  where  /j,q  denotes  the  magnetic  perme¬ 
ability,  and  form  the  Gibbs  energy  relations 

G(H,  M,  T)  =  ip(M,  T)  -  hqHM  (26) 

or 

G{H ,  M,  T)  =  i](M,  T)  -  HM  (27) 

by  incorporating  hq  int°  V’-  For  increasing  H ,  the  behavior  of  G  with  if  given  by  (25)  is  depicted 
in  Figure  4(a).  In  the  absence  of  anisotropic  effects  or  applied  stresses,  G  approximates  the  energy 
landscape  exhibited  at  the  lattice  level  in  homogeneous  materials. 


3  Local  Average  and  Anhysteretic  Magnetizations 

3.1  Local  Magnetization 

For  conditions  in  which  thermal  after-effects  [2]  are  negligible,  the  local  average  magnetization  M 
at  the  lattice  level  is  determined  by  minimizing  the  Gibbs  relations  (26)  or  (27)  whereas  the  Gibbs 
energy  must  be  balanced  with  the  thermal  energy  through  Boltzmann  principles  if  thermal  effects 
are  significant.  We  consider  these  two  regimes  in  Sections  3.1.1  and  3.1.2  and  then  illustrate  in 
Section  3.1.3  that  the  model  which  incorporates  thermal  energy  limits  to  the  case  of  no  thermal 
activation  when  reference  volumes  V  are  taken  to  be  arbitrarily  large. 


3.1.1  Negligible  Thermal  Effects 

For  conditions  in  which  thermal  after-effects  are  negligible,  the  local  average  magnetization  M  is 
determined  from  the  necessary  conditions 


dG  _  d2G 
i).\I  ~  ’  <).\I- 


(28) 


When  the  statistical  mechanics  relation  (18)  is  employed  for  the  Helmholtz  energy,  this  yields  the 
Ising  relation 

M(H)  =  Ms tanh  (29) 

where  a  and  a(T )  are  defined  in  (22).  The  behavior  of  the  kernel  or  hysteron  is  illustrated  in 
Figure  5(a). 


Remark  1. 


The  Ising  relation  (29),  whose  input  is  the  effective  field 


He  =  H  +  aM, 


(30) 


is  fundamental  in  a  number  of  hysteresis  models  for  ferromagnetic  and  ferroelectric  compounds. 
This  relation  was  directly  employed  for  quantifying  the  anhysteretic  component  of  unified  models 
developed  in  [54],  Furthermore,  it  is  illustrated  in  [54,  55]  that  if  one  relaxes  the  constraint  that 
moments  have  only  the  orientations  a  *  =  ±1  and  considers  uniformly  distributed  moments,  one 
obtains  the  Langevin  relation  M  =  C{He)  =  Ms[coth(He/a)  —  a/He ]  which  agrees  with  the  Ising 
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(a)  (b) 

Figure  5:  (a)  Limiting  kernel  (29)  obtained  with  the  Helmholtz  energy  (18)  in  the  absence  of  thermal 
activation,  (b)  Limiting  kernel  (32)  provided  by  the  Helmholtz  relation  (25). 


relation  M  =  Mstanh(He/a)  through  first  order  terms  -  see  [2,4]  for  a  derivation  of  the  Langevin 
equation  in  the  context  of  magnetic  materials.  The  Langevin  model  M  =  MsC(He),  with  He  specified 
by  (30),  is  employed  when  quantifying  the  anhysteretic  magnetization  in  the  domain  wall  theory  of 
Jiles  and  Atherton  [10,  11]  as  well  as  the  transducer  models  based  on  that  theory  [45-47]  —  e.g., 
see  (1).  Finally,  translates  of  the  Ising  relation  r(x)  =  tanh(x)  form  appropriate  ridge  functions 
for  generalized  Preisach,  or  KrasnoseTskn-Pokrovskii  characterizations  [56-58].  Hence  this  relation 
plays  a  fundamental  role  in  two  of  the  present  macroscopic  theories  outlined  in  Section  1. 

The  local  average  magnetization  M  resulting  from  (25)  is  elementary  in  the  sense  that  it  is 
piecewise  linear  but  is  complicated  by  the  fact  that  a  history  of  moment  switches  must  be  maintained 
to  ascertain  which  branch  of  the  hysteron  is  active.  Enforcement  of  the  necessary  condition(28)  yields 

M=1H  +  Mr6  (31) 

T] 

where  <5=1  for  positively  oriented  moments  and  5  =  —  1  for  negative  orientations.  To  quantify  6 
in  terms  of  initial  moment  configurations  and  previous  switches,  we  employ  Preisach  notation  - 
e.g.,  see  [33,56,58]  —  and  take 


[M(H;Hc,Om=  < 


(  [M(H;Hc,O}(0)  ,r(t)  =  0 

§  -  Mr  ,  r(t)  0  and  H (max  r(t))  =  -Hc 

,  r(t)  0  and  H (max  r(t))  =  Hc. 


I  f  +  Mr 


Here 


[M(H;Hc,Om  =  < 


H 


H 


-  Mr  ,  H( 0)  <  -Hc 

,  -Hc  <  H( 0)  <  Hc 
+  Mr  ,  H( 0)  >  Hc 


denotes  the  initial  moment  distribution  and  transition  times  are  designated  by 

r(t)  =  {t  G  (0,  tf]  |  H(t)  =  —Hc  or  H(t)  =  Hc} 

where  tf  denotes  the  final  time  under  consideration. 

The  dependence  of  M  on  the  local  coercive  field 


(32) 


(33) 


(34) 


Hr  =  ?7 (Mr  -  Afj) 


(35) 


is  indicated  as  a  prelude  to  the  discussion  in  Section  4  where  Hc  is  assumed  distributed  to  accom¬ 
modate  material  nonhomogeneities. 
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The  behavior  of  the  limiting  kernel  (31)  or  (32)  is  compared  with  its  statistical  mechanics  coun¬ 
terpart  (29)  in  Figure  5.  It  is  observed  that  the  primary  difference  between  the  kernels  occurs  in 
the  saturation  behavior  at  high  fields.  The  kernel  (32)  predicts  a  linear  relation  between  H  and 
M  after  moment  switching  whereas  the  kernel  (29)  exhibits  saturation  behavior  to  a  local  average 
magnetization  value  Ms . 

Remark  2.  A  comparison  of  the  upper  and  lower  branches  of  the  Ising  kernel  plotted  in  Figure  5(a) 
illustrates  that  this  kernel,  obtained  from  the  energy  relation  (18),  yields  noncongruent  behavior. 
Hence  it  can  be  used  to  characterize  the  noncongruency  measured  in  certain  operating  regimes. 

Remark  3.  The  limiting  kernel  (32)  provides  reversible  behavior  at  high  fields  due  to  the  fact  that 
the  kernel  does  not  saturate. 


3.1.2  Thermal  After-Effects 

To  incorporate  the  thermal  mechanisms  which  produce  phenomena  such  as  after-effects  [2],  it  is 
necessary  to  balance  the  Gibbs  energy  G  with  the  relative  thermal  energy  kT/V,  over  the  reference 
volume  V ,  through  the  Boltzmann  density  relation 

p{G)  =  Ce~GV/kT  (36) 

which  specifies  the  probability  of  attaining  an  energy  level  G  for  a  fixed  field  input.  The  constant  C  is 
chosen  to  ensure  a  probability  of  unity  when  p  is  integrated  over  all  admissible  moment  configurations. 


Gaussian  Behavior  of  p 

To  illustrate  the  behavior  of  p  for  the  Gibbs  energy  G  =  ip  —  HM  constructed  using  the  piecewise 
Helmholtz  model  (25),  we  consider  the  specific  energy  profile  depicted  in  Figure  6  for  which  it  is 
assumed  that  H  >  0  and  G(M+in)  <  G{M(nin )  <  G(Mq).  The  relative  minima 


M-JH)  =  |  -  Mr, 

K,nW)  =  f  + 


(37) 


result  from  the  necessary  condition  (28)  utilized  when  constructing  the  limiting  model  (31)  or  (32). 
The  local  coercive  field  Hc  for  which  M~-  =  —Mj  =  Mq  is  given  by  (35). 

From  (25),  it  follows  that  for  M  <  —  Mj,  the  Boltzmann  probability  can  be  formulated  as 


p(G(H,M )) 


C(T)e~G(H’M)v/kT 

e-[§»7  {M+MR)2-HM]V/kT 


r — Mj 


-[±v(M+MR)2-HM}V/kTdM 


e-(M-M-in)2vV/2kT 


-M, 


o-{M-M-in)^V/2kTdM 


C(T,0)e 


—(M—M~ .  )2/2/32 

v  min '  '  ^ 


(38) 
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Figure  6:  (a)  Gibbs  energy  profile  with  a  high  level  ( - )  and  low  level  ( - )  of  thermal  activation 

in  the  Boltzmann  probability  =  Ce~GV^kT .  (b)  Local  magnetization  M  given  by  equation  (49) 

with  high  thermal  activation  ( - )  and  limiting  magnetization  M  specified  by  (32)  in  the  absence 

of  thermal  activation  ( - ). 


where 


r»  — Mj 


o-{M-M~in)2W2{T)dM 


C(T,f5)  = 

Similarly,  for  M  >  M j ,  the  probability  has  the  form 

e-(M-MLn)2W/2  kT 

M(M)  =  — 


-l 


Mj 


e-(M-M+in)^V/2kTdM 


(39) 


(40) 


Relations  (38)  and  (40)  illustrate  the  Gaussian  behavior  of  the  Boltzmann  probabilities  for  the 
piecewise  quadratic  Helmholtz  function  if)  while  (39)  illustrates  that  the  variance  (5 2  is  proportional  to 
the  relative  thermal  energy  kT/V.  From  a  physical  perspective,  low  relative  thermal  energy  implies 
that  fewer  moments  achieve  the  energy  required  to  overcome  energy  barriers  thus  producing  steep 
transitions  in  the  local  relation  between  H  and  M . 

To  illustrate  the  Dirac  nature  of  n(G)  in  (38)  as  kT/V  decreases,  let  j  =  1/(5  and  define  the 
sequence 


C{T,j)e-{M~Mmin?f/2  ,  M  <  —  Mj 

0  ,  M  >  -Mj  . 


(41) 


The  sequence  {4>j}  satisfies  properties  (i)-(iii)  of  Theorem  1  in  7  and  hence  constitutes  a  Dirac  family. 
It  follows  immediately  that 


lim  u(M )  =  lim  6AM) 

kT/v^O  j- oo 


(42) 


=  5(M-M~J. 


Analogous  behavior  is  exhibited  at  M+m  as  depicted  in  Figure  6(b). 


Transition  Likelihoods  and  Local  Average  Magnetization 

Because  the  Boltzmann  relation  (36)  quantifies  the  balance  between  the  Gibbs  and  relative  ther¬ 
mal  energies,  it  is  employed  when  modeling  the  fraction  of  positively  and  negatively  oriented  mo¬ 
ments,  the  average  magnetizations  due  to  the  two  configurations,  and  the  likelihoods  that  moments 
changes  configurations  for  a  given  input  field  level. 
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Recalling  that  N _  and  N+  respectively  denote  the  number  of  negatively  and  positively  oriented 
moments,  we  denote  the  respective  moment  fractions  by  x_  =  N_/N  and  x+  =  N+/N  where 


X-  +  x+  =  1  (43) 

since  N—  +  N+  =  1.  The  evolution  of  moment  fractions  is  determined  by  the  differential  equations 


X+  =  ~P- 1 _ x_|_  +p _ |_x_ 

X-  =  —p _ yX-  +  P- 1 _ X+ 


(44) 


which  can  be  simplified  to 

x+ = -p+-X+ +  p-+(  1  —  x+)  (45) 

through  the  identity  (43).  For  a  demagnetized  material,  initial  conditions  can  be  taken  to  be  x_  = 
x+  =  l/2. 

The  expected  or  average  magnetizations  due  to  negatively  and  positively  oriented  moments  are 
defined  by 

/•Mo(T)  poo 

{MJ)  =  /  Mp(M)  dM  ,  (M+)  =  /  Mp(M)  dM  (46) 

4-oo  Jm0{T) 

where  Mq{T )  denotes  the  unstable  equilibrium  of  G  as  depicted  in  Figure  6(a).  For  the  piecewise 
quadratic  Helmholtz  energy  functional  (25),  the  evaluation  of  the  integrals  in  (46)  is  simplified  by 
replacing  the  limit  Mq(T )  respectively  by  —Mj  and  M j  in  the  definitions  of  {MJ)  and  (M+) .  This  can 
be  motivated  by  observing  that  maximum  restoring  forces  occur  at  the  inflection  points  as  detailed 
on  pages  332-333  of  [3]  or  pages  486-487  of  [2].  Furthermore,  these  points  coincide  in  the  limit  of 
negligible  thermal  activation  as  illustrated  in  Section  3.1.3.  With  this  approximation,  we  have 


f*—Mj 


Me-G{H,M)V/kTdM 


Me~G(H,M,T)V/kT  dM 


(M-)  = 


e~G(H,M)V/kT  dM 


(M+)  = 


lM j 


/  e-G(H,M,T)V/kTdM 

lM ! 


(47) 


The  likelihood  of  switching  from  a  positive  moment  orientation  to  negative,  and  conversely,  are 
respectively  quantified  by 


P+-  = 


r*Mj 


'  Mj — € 


-G(E,M)V/kT  dM 


T(T) 


/ 
J  A 


-G(E,M)V/kT  dM 


P-+  = 


r—  Mj+e 


'  —Mj 


-G(E,M)V/kT  dM 


TIT )  r~Mi+e 

{  ’  e~G^v/kTdM 


(48) 


Mj-e 


where  e  is  taken  to  be  a  small  positive  constant.  The  quotient  of  integrals  is  a  probability  and  hence 
is  unitless.  The  relaxation  time  T  is  the  reciprocal  of  the  frequency  at  which  moments  attempt  to 
switch  so  y  has  units  of  ^ .  This  yields  the  correct  units  in  the  differential  equations  (44)  and 
(45).  Moreover,  we  note  that  T 2  is  considered  to  be  inversely  proportional  to  the  relative  thermal 
energy  so  that  T (T)  =  T\  sj V/kT ;  hence  increased  temperature  lead  to  increased  thermal  relaxation 
behavior.  For  materials  having  a  single  relaxation  time,  T\  is  constant  whereas  variable  relaxation 
times  may  need  to  be  identified  for  materials  exhibiting  distributed  relaxation  behavior. 

With  the  moment  fractions,  expected  magnetization  values,  and  transition  likelihoods  thus  de¬ 
fined,  the  local  average  magnetization  for  the  lattice  is 


M  =  x+  (M. |_)  +  x-  {MJ)  . 


(49) 
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The  behavior  of  the  local  model  (49),  which  incorporates  thermal  after-effects  (thermal  relaxation) 
is  compared  in  Figure  6(b)  with  the  relation  (32)  obtained  by  simply  minimizing  the  Gibbs  energy. 
For  values  of  kT /V  on  the  order  of  G,  a  significant  number  of  moments  achieve  the  relative  thermal 
energy  required  for  switching  in  advance  of  the  local  coercive  field  Hc.  This  produces  a  smooth 
transition  between  the  limiting  minima  (37)  of  the  hysteron.  For  diminishing  values  of  kT /V  as 
compared  with  G,  fewer  moments  achive  the  thermal  energy  required  for  pre-coercive  switching 
which  produces  increasingly  steep  transitions  between  orientations.  This  convergence  is  rigorously 
established  in  the  next  section. 


3.1.3  Limiting  Behavior  of  Local  Magnetization  Model 

The  local  model  (49)  incorporates  thermal  after-effects  by  employing  the  Boltzmann  relation  (36) 
to  balance  the  Gibbs  energy  G  and  relative  thermal  energy  kT /V  whereas  the  local  model  (32)  was 
derived  in  the  absence  of  thermal  relaxation  mechanisms  simply  by  minimizing  the  Gibbs  energy. 
We  rigorously  establish  here  the  convergence  of  (49)  to  (32)  in  the  limit  kT /V  — ►  0  of  increasing 
control  volumes  and  hence  diminishing  relative  thermal  energies.  To  clarify  the  discussion,  we  con¬ 
sider  the  representative  energy  landscape  depicted  in  Figure  6(a)  -  however,  the  analysis  techniques 
accommodate  general  energy  configurations. 

We  consider  first  the  convergence  of  the  expected  magnetization  relations  (47).  For  negative 
moments,  we  consider  the  Dirac  sequence  {4>j}  defined  in  (41),  define  the  function  f(M)  =  M,  and 
consider  the  interval  [a,  b]  =  [M~m,  My]  where  Mmm  is  defined  in  (37).  Hence  /  is  continuous  on  M 
and  satisfies  the  decay  property  (iv)  in  Theorem  1  of  7.  It  then  follows  that 


lim  (M_) 

kT/V^O 


lim  /  M(/>j(M  -  M-JdM 

3-+°°J_0O 


m: 


(50) 


Analogous  arguments  can  be  used  to  demonstrate  that  (M+)  —>  M+in  as  kT /V  — >  0. 

To  illustrate  the  convergence  of  the  transition  likelihoods  for  a  fixed,  relaxation  time  T  (T) ,  we 
modify  the  sequence  {4>j}  defined  in  (41)  for  the  interval  (— oo,  —Mj  +  e].  The  function  /  is  specified 
to  be 


f(M) 


0  ,  M  <  -Mj 
1  ,  M  >  -Mj 


(51) 


and  the  interval  [a,  b]  is  taken  to  be  [— 2Mmin,  —  Mj]  or  [—Mi,  Mo].  Since  /  again  satisfies  (iv)  in 
Theorem  1  of  7,  we  obtain  the  convergence 


lim  p _ l 

kT/V^O 


lim  1  [°°  /(M)^-(M-  M-JdM 

J^CX)  T  (1  )  J_OQ 


H  <HC 
H  >  Hc 


(52) 


for  Hc  defined  by  (35).  Similar  analysis  for  positively  oriented  moments  in  the  considered  energy 
landscape  yields 

lim  p+_  =  0.  (53) 

kT/V^O 

We  let  C+  =  x+(0)  and  C-  =  x_(0),  (+  +  C-  =  1,  denote  the  initial  moment  fractions  and  consider 
the  behavior  of  the  differential  equation  (45)  governing  the  evolution  of  x+.  Under  the  assumption 
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that  H,  which  is  parameterized  with  respect  to  time,  is  increasing,  we  let  t  =  tc  denote  the  time  at 
which  H(t)  =  Hc.  The  solution  to  (45)  in  the  limit  kT/V  — >  0  is  then 


x+(t)  = 


C+ 


t  <  tr 


[M(H)\(t)  = 


[M(H)](t)  = 


[  1  -  (1  -  C+)e-(*-tc)/r 

,  t  >  tc 

ion  is 

^  +  (2C+  -  1  )MR 

,  t  <tc 

MI)  +  [1  _  2(1  -  C+)e-(^ 

C)/T]  Mr  ,  t  >  tc 

limits  to 

1  m  +  (2(+-i)MR  , 

t  <  tc 

1  ^  +  MR 

t>tc 

c~[M-in(Hm + c+[M+in(Hm 

iKinimt) 


t  <  tc 
t  >  tc 


which  is  precisely  (32). 


(54) 


(55) 


(56) 


Remark  4.  To  summarize,  the  linear  kernel  (31)  or  (32)  can  be  accurately  employed  when  thermal 
effects  are  negligible  (kT/V  is  small)  and  relaxation  times  T  are  small  compared  with  drive  frequen¬ 
cies.  Otherwise,  one  should  employ  the  kernel  (49)  or  an  asymptotic  relation  of  the  form  (55)  to 
accommodate  thermal  activation  or  long  relaxation  times. 


3.2  Local  Anhysteretic  Magnetization 

The  relations  (32)  and  (49)  characterize  the  local  hysteretic  H-M  behavior  at  the  lattice  level  when 
the  piecewise  quadratic  relation  (25)  is  used  to  quantify  the  Helmholtz  energy.  The  energy  framework 
used  to  establish  these  relations  also  quantifies  the  anhysteretic  H-M  behavior  which  is  experimen¬ 
tally  achieved  by  applying  sufficiently  large  AC  fields  superimposed  on  a  DC  bias  field.  From  a 
theoretical  perspective,  the  local  anhysteretic  magnetization  Man  represents  the  locus  of  magnetiza¬ 
tion  values  which  would  occur  in  materials  devoid  of  inclusions.  It  can  also  be  theoretically  formulated 
as  the  magnetization  achieved  when  relaxation  times  T (T)  are  sufficiently  small  compared  with  drive 
frequencies  that  moments  achieve  global  equilibria. 

We  illustrate  the  latter  theoretical  interpration  in  the  context  of  the  local  magnetization  model 
(49)  derived  under  the  assumption  that  G  and  kT/V  are  balanced  through  the  relation  (36).  The 
condition  of  moment  equilibrium  yields  x+  =  x_  =  0  in  (44)  which  in  turn  implies  that  equilibrium 
solutions  and  x _  satisfy  the  relation 


To  demonstrate  the  implication  of  (57) 

(48),  it  follows  immediately  that  p _ |_  =  p _| _ and  hence  x+  =  x_.  From  the  conservation  relation 

x+  +  X-  =  1,  it  is  deduced  that  x+  =  X-  =  2 ,  regardless  of  the  initial  conditions  x+(0)  and  x_(0). 
The  rate  at  which  the  relations  converge  to  equilibrium  values  is  determined  by  the  relaxation  time 
T(T),  with  smaller  values  of  T  producing  more  rapid  equilibration. 


x± 

X- 


P-+ 

P-+' 


(57) 


consider  first  the  case  when  H  =  0.  From  the  definition 
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To  ascertain  the  resulting  anhysteretic  magnetization  for  H  =  0,  we  note  that  the  symmetry  of 
(47)  implies  that 

(M+)  =  -(M_)  (58) 

at  equilibrium.  When  combined  with  the  fact  that  =  x_  =  ,  (49)  yields 

Man{H  =  0)  =  0.  (59) 

The  field-dependence  of  p.\ _ ,p _ (M+)  and  (M_)  precludes  a  similar  exploitation  of  symmetries 

for  H  0.  However,  Man(H )  can  be  easily  computed  by  numerically  approximating  (49)  with 
sufficiently  small  T  -  recall  that  lo  =  j*  quantifies  the  frequency  at  which  moments  attempt  to 
switch.  The  Gibbs  energy  G  at  the  field  value  Hq  =  2000  A/m,  unnormalized  density  p{G)  = 
e-GV/kT ^  ancj  resuitmg  local  anhysteretic  magnetization  obtained  with  T  =  1.0  x  10-13  sec,  and 
relative  thermal  energies  kT  /V  =  5.0  x  105  and  kT /V  =  7.14  x  106  are  plotted  in  Figure  7. 

It  is  observed  that  when  kT /V  is  significant  compared  with  G,  thermal  fluctuations  produce 
switching  between  wells  thus  yielding  a  gradual  anhysteretic  transition  between  positive  and  negative 
saturation  magnetizations.  As  kT /V  becomes  increasingly  small,  the  local  anhysteretic  magnetiza¬ 
tion  M an  provided  by  (49)  converges  to 


M(H)  =  —  +  Mr5 

V  (60) 

(5  =  sign  (H). 

The  limiting  relation  (60)  can  be  interpreted  as  the  locus  of  magnetization  values  which  would  occur 
in  the  absence  of  inclusions  —  which  is  manifested  by  Hc  =  0  in  the  local  model. 

Remark  5.  We  note  that  the  local  magnetization  relation  (29)  obtained  from  the  necessary  condition 
jnfj  =  0,  with  the  statistical  mechanics  model  for  G,  also  yields  an  anhysteretic  magnetization  for 


Legend 

- kT/V=  5.0  x10s 

- kT/V  =  7.1x106 


Figure  7:  (a)  Gibbs  energy  G,  and  (b)  unnormalized  Gaussian  densities  l-i(G)  =  e  GV/kT  with 
Hq  =  2000  A/m.  (c)  Anhysteretic  magnetization  Man  given  by  (49)  with  T  =  1.0  x  10”13  sec. 
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certain  values  of  a  and  a.  This  is  analogous  to  the  Jiles-Atherton  framework  in  which  the  Langevin 
expression  (1)  is  employed  to  quantify  Man.  Due  to  its  simplicity,  however,  we  generally  employ 
(60)  as  a  kernel  when  characterizing  the  anhysteretic  magnetization  for  regimes  in  which  thermal 
relaxation  is  negligible. 

4  Macroscopic  Magnetization  Models 

The  local  magnetization  models  (29,  (32)  or  (49)  were  derived  by  constructing  appropriate  energy 
relations  at  the  lattice  level.  For  homogeneous  materials  with  uniform  effective  fields,  these  relations 
hold  throughout  the  material  and  hence  will  also  provide  macroscopic  models.  By  construction,  they 
exhibit  the  steep  transitions  depicted  in  Figures  4,  5  and  6  since  moments  are  assumed  to  switch 
instantaneously  once  they  achieve  the  energy  required  to  overcome  energy  barriers.  The  models 
in  this  form  prove  adequate  for  characterizing  the  hysteretic  behavior  of  certain  materials  which 
exhibit  small  demagnetizing  factors,  certain  uniaxial  wires  and  films  or  annealed  toroidal  specimens. 
Hysteresis  loops  exhibiting  nearly  instantaneous  transitions  due  to  these  factors  are  illustrated  for 
a  uniaxial  nickel-iron  film,  a  magnetically  annealed  core  of  cobalt  ferrous  ferrite  and  manganese- 
magnesium  ferrite  on  page  298  of  Craig  and  Tebble  [59].  However,  the  transitions  provided  by  these 
local  models  are  too  steep  to  provide  accurate  characterization  of  general  polycrystalline  magnetic 
materials.  To  extend  the  local  models,  we  consider  certain  parameters  in  the  models  to  be  statistically 
distributed  to  reflect  variations  in  the  lattice  structure,  exchange  energies  and  grain  orientations.  The 
resulting  macroscopic  magnetization  models  accurately  characterize  both  major  and  biased  minor 
loops  in  a  wide  range  of  ferromagnetic  materials. 

4.1  Statistical  Mechanics  Model 

An  implicit  assumption  made  when  deriving  the  Helmholtz  energy  relation  (18)  used  to  construct  the 
local  magnetization  models  (29)  and  (49)  is  that  the  exchange  energy  4>o  is  constant  throughout  the 
lattice.  This  implies  that  the  bias  field  Hh  =  2vm  and  Curie  temperature  Tc  =  are  constant  which 
yields  a  constant  mean  field  coefficient  a  =  and  constant  coefficient  a(T )  =  in  the  models 
(29)  and  (49).  However,  for  nonhomogeneous,  polycrystalline  materials  with  variable  magnetization, 
this  assumption  is  overly  simplistic.  For  such  materials,  it  is  more  reasonable  to  assume  instead  that 
$0  is  statistically  distributed  which  motivates  the  consideration  of  statistically  distributed  parameters 
in  the  macroscopic  magnetization  models.  Additionally,  material  nonhomogeneities,  variable  grain 
orientations,  nonuniform  stress  distributions,  and  variations  due  to  texture  motivate  consideration 
of  statistically  distributed  model  parameters. 

Because  <f>o  quantifies  the  energy  required  to  reorient  a  moment  when  the  lattice  is  completely 
ordered,  the  assumption  that  $o  is  statistically  distributed  implies  that  the  exchange  energy  between 
spins  or  moments  is  distributed.  Through  (12),  this  implies  that  the  exchange  integral  J  is  variable 
rather  than  constant  as  assumed  for  homogeneous  materials.  At  the  quantum  level,  the  variability 
of  Jij  employed  in  (5)  is  incorporated  by  modeling  the  overlap  of  electron  wave  functions  whereas 
at  the  macroscopic  level,  it  is  incorporated  by  considering  Hh  and  a  to  be  statistically  distributed. 
We  consider  the  construction  of  macroscopic  mean  field  models  which  accommodate  the  microscopic 
variations  in  the  exchange  integral  and  lattice  energy  4>o- 

We  first  make  the  assumption  that  $0  is  normally  distributed  about  a  mean  value  of  4>o-  If  hi,  V 
and  Ms  remain  constant,  the  bias  field  Hh  =  will  then  be  normally  distributed  with  mean  Hh- 
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A  resulting  macroscopic  magnetization  model  is 

/OO  _  2 

[M(H;Hh,m)e-{Hh-Hh>  ^  dHh  (61) 

-OO 

where  b  and  C  are  constants  and  M  is  specified  by  (29)  or  (49).  Because  a  =  this  is  equivalent 
to  employing  effective  fields 

He  =  H  +  aM  (62) 

where  a  is  normally  distributed.  This  should  be  compared  with  a  number  of  current  hysteresis 

models  which  employ  effective  fields  of  the  form  (62)  with  fixed  a  (e.g.,  [4, 11,22]). 

More  general  effective  fields  can  be  incorporated  if  it  is  assumed  that  in  addition  to  variations 
in  Hh,  the  effective  field  He  at  the  domain  level  is  normally  distributed  about  the  applied  field  H. 
Since  He  =  H  +  Hj,  where  Hi  denotes  the  interaction  field,  this  yields  the  macroscopic  model 

/oo  roo  _  2 

/  \M(H  +  Hi;Hh,li)]{t)e-H^2b2e-(Hh-Hh)  /2b2dHKiHh  (63) 

-oo  J  —  OO 

where  b 2  determines  the  variance  about  H. 


4.2  Piecewise  Quadratic  Helmholtz  Relations 


For  the  model  derived  using  the  piecewise  quadratic  Helmholtz  energy  (25),  we  incorporate  lattice 
variations  due  to  nonhomogeneities,  polycrystallinity  and  material  inclusions  by  considering  the  local 
coercive  field  Hc  to  be  stochastically  distributed.  To  enforce  Hc  >  0,  we  make  the  initial  assumption 
that  Hc  satisfies  a  lognormal  distribution  with  density 


vi(Hc)  =  ci  exp 


In  {Hc/Hr2' 
2c 


(64) 


where  c,  ci  and  Hc  are  positive  constants.  It  is  illustrated  in  [9],  where  this  distribution  is  considered 
in  the  context  of  Preisach  models,  that  if  c  is  small  compared  with  Hc ,  the  mean  and  variance  of  / 
have  the  approximate  values 

(Hc)  «  Hc  ,  a^2Hcc.  (65) 

The  macroscopic  magnetization  model  based  on  this  distribution  of  Hc  is 

POO 

[M(H)}(t)=  /  \M(H;Hc,0}(t)vi(Hc)dHc  (66) 

Jo 

where  M  is  given  by  (32)  or  (49).  The  relations  (65)  can  be  used  to  obtain  initial  parameter 
estimates  from  attributes  of  measured  data.  For  certain  materials,  the  distribution  of  Hc  can  be 
taken  as  Gaussian;  however,  positivity  should  still  be  enforced  in  the  integration  limits. 

To  incorporate  variability  in  the  exchange  energy  $o>  we  consider  the  effective  field  He  to  be  nor¬ 
mally  distributed  about  the  applied  field  H.  When  combined  with  (66),  this  yields  the  macroscopic 
magnetization  model 

/*  oo  /*oo  _ 

[M(H)](t)  =  C  /  \M(H  +  HI-Hc,C)]{t)e-HV2h\-^H^HJ^dHIdHc  (67) 
J  0  J — oo 

where  C  and  b  are  positive  constants  and  M  is  given  by  (32)  or  (49).  When  the  kernel  M  is  computed 
using  (49),  the  model  incorporates  certain  relaxation  mechanisms  including  magnetic  after-effects. 
However,  it  does  not  incorporate  elastic  effects  or  eddy  current  dynamics  in  this  formulation  so  it 
should  be  employed  in  low  frequency  regimes. 
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4.3  General  Density  Formulation 

The  macroscopic  model  formulations  (63)  and  (67)  are  based  on  a  priori  assumptions  regarding 
the  normal  or  lognormal  nature  of  underlying  densities.  In  certain  cases,  the  normal  behavior  of 
parameters  can  be  argued  using  statistical  theory  based  on  the  central  limit  theorem ;  e.g.,  see  [9]. 
In  general,  however,  the  choice  of  normal  or  lognormal  distributions  is  based  solely  on  mathematical 
attributes  rather  than  physical  or  energy  principles.  These  mathematical  assumptions  can  be  avoided 
by  formulating  the  macroscopic  models  in  terms  of  general  densities  to  be  estimated  through  a  least 
squares  fit  to  data. 

To  illustrate,  let  f\  and  ^2  designate  the  densities  respectively  associated  with  local  coercive  and 
effective  fields.  To  satisfy  physical  criteria,  we  assume  that  v \  and  z^2  satisfy  the  conditions 


(i)  v\(x)  defined  for  x  >  0, 

(ii)  V2(~x)  =  v2(x), 

(Hi)  Ms)  |  <  Cle~aix  ,  Ms)  |  <  C2e~a2X 


(68) 


for  positive  ai,ct2,ci,C2.  The  restricted  domain  in  (i)  reflects  the  fact  that  the  coercive  field  Hc 
is  positive  whereas  the  symmetry  enforced  in  the  effective  field  through  (ii)  yields  the  symmetry 
observed  in  low- field  Rayleigh  loops.  Hypothesis  (iii)  incorporates  the  physical  observation  that  the 
coercive  and  interaction  fields  decay  as  a  function  of  distance  and  guarantees  that  integration  against 
the  piecewise  linear  kernel  yields  finite  magnetization  values. 

For  the  piecewise  quadratic  Helmholtz  energy  (25),  The  resulting  macroscopic  magnetization 
model  is  then  given  by 


roo  roc 

[M{H)](t)  =  /  /  v\(Hc)v2(Hi)\M{H  +  Hi;  Hc,£)](t)  dHi  dHc 

J  0  J —00 


(69) 


V(HC,  Ht)[M(H  +  Hr,  Hc,  £)](£)  dH:  dHc 


'  0  J  —00 


where  M  is  specified  by  (31)  or  (49).  Formulation  in  terms  of  the  product  density  v  is  more  general 
whereas  retention  of  the  components  v\  and  v2  can  facilitate  subsequent  implementation. 

Remark  6.  A  comparison  of  (69)  with  (2)  indicates  the  manner  through  which  the  framework 
provides  an  energy  basis  for  certain  extended  Preisach  models  as  detailed  in  [32,  33].  Techniques 
analogous  to  those  developed  for  Preisach  models  [27]  have  been  developed  to  identify  the  general 
density  v  in  an  efficient  manner  [60,  61],  thus  exploiting  the  similarity  between  the  two  frameworks. 
As  detailed  in  Sections  1  and  5,  the  formulation  (69)  is  advantageous  over  classical  Preisach  models 
in  the  manner  through  which  it  relaxes  reversibility,  deletion,  and  congruency  criteria  and  incorpo¬ 
rates  temperature  and  rate- dependencies  in  the  basis  M  rather  than  in  the  parameters  u.  For  these 
reasons,  it  provides  an  energy  basis  for  certain  extended  Preisach  models  [9]. 


4.4  Anhysteretic  Magnetization  Model 

It  was  illustrated  in  Section  3.2  that  for  the  piecewise  quadratic  Gibbs  energy,  the  local  anhysteretic 
magnetization  Man  followed  naturally  from  (49)  for  thermally  active  regimes  or  (60)  in  the  absence 
of  thermal  after-effects  or  relaxation.  The  global  anhysteretic  model  follows  directly  from  the  general 
hysteresis  model  (67)  or  (69)  but  can  be  simplified  substantially  since  coercive  fields  play  no  role  in 
the  anhysteretic  material  behavior. 
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We  consider  first  the  anhysteretic  model  for  the  a  priori  density  choices 


=  ?Le-MHc/Hc)/2c[* 

h 

ru  \  —  ^2  „-H2T/2b2 

',2{H,)-wse 

where  c\  =  c\/I\  and  C2  =  C2/b\f2n  are  expressed  in  terms  of  the  normalization  constants 


h  = 


roo  _  roo 

/  e-MHc/Hc)M 2  dHc  ,  bV^  =  /  e~HV2b"  dHj. 

J  0  J — oo 


The  anhysteretic  magnetization  is  then  given  by 


~  ~  roc  roc  _ 

Man{H)  =  -±.—?=  /  M(H  +  HI)e-[HHc/Hc)/2c\2e~Hl/2b2  dHrdHc 

h  bV 2n  Jo  J- oo 


h  b\/2iT 

/OO 

M{H  +  Hj)e-^ 

-OO 


dHf 


where 


C  = 


cic2 

6v/2vr 


(70) 


(71) 


(72) 


(73) 


If  solely  quantifying  anhysteretic  material  behavior,  one  can  treat  the  constant  C  as  a  material 
parameter  to  be  identified  whereas  if  correlating  modeled  hysteretic  and  anhysteretic  properties,  one 
should  identify  the  constants  c\  and  C2.  In  the  absence  of  thermal  after-effects,  the  local  relation 
(60)  yields 

M(H  +  77/)  =  —  +Hl  +  Mr6 

d  (74) 

S  =  sign  (H  +  Hi) 

whereas  the  kernel  (49)  can  be  employed  if  thermal  activation  is  significant. 

Additional  generality  can  be  obtained  through  the  formulation 


Man(H)  =  /  M(H  +  Hj)v2(Hj)  dHT 


(75) 


where  v-2  is  a  general  density  satisfying  the  assumption  (68).  As  with  the  parameterized  formulation, 
normalization  constants  must  be  accounted  for  if  comparing  the  anhysteretic  model  (75)  and  hys¬ 
teresis  model  (69)  —  this  reflects  the  price  paid  for  employing  unnormalized  density  formulations  to 
simplify  notation. 


4.5  Model  Implementation 

Two  issues  must  be  addressed  when  implementing  the  hysteresis  models  (63)  or  (67);  (i)  approxima¬ 
tion  of  the  integrals,  and  (ii)  efficient  implementation  of  the  conditional  relations  (32).  Because  both 
issues  are  crucial  for  providing  algorithms  that  permit  efficient  system  design  and  real-time  control 
implementation,  we  summarize  pertinent  details.  For  simplicity,  we  focus  on  the  implementation  of 
(67)  and  note  that  analogous  constructs  exist  for  (63). 
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4.5.1  Quadrature  Techniques 

The  integrals  can  be  approximated  either  on  the  original  infinite  and  semi-infinite  domains  or  on  finite 
domains  determined  by  the  exponential  decay  properties  of  the  integrands.  On  the  infinite  domain, 
Gauss-Hernrite  quadrature  formulae  apply  whereas  Gauss-Laguerre  points  and  weights  apply  for 
the  semi- infinite  integrals  [62].  As  illustrated  in  Figures  8(a)  and  8(c),  the  exponential  decay  of 
the  densities  can  also  be  employed  to  determine  finite  intervals  where  Gauss-Legendre  formulae  are 
accurate. 

In  all  cases,  approximation  of  (67)  yields 


Ni  Nj 

\  \  ^ 


[M(H)}(t)  =  C^^[M(HI.+H-,HCi,Zj)}(t)e 


—Hj  /2b2 


[\n(HcJHc)/2c\ 


ViWj 


(76) 


i= 1  j= 1 


where  Hjj ,  HCi  denote  the  abscissas  and  Vi ,  Wj  are  weights  associated  with  the  respective  quadrature 
formulae.  At  H  =  0,  M  =  0,  the  initial  moment  distribution  corresponds  with  the  quadrature 
points  as  illustrated  in  Figure  8(b). 

To  further  illustrate,  we  consider  the  construction  of  Gauss-Legendre  points  and  weights  on  the 
interval  [— L,L\  using  a  4  point  composite  quadrature  rule.  On  each  subinterval  [hj-i,hj\,  where 
hj  =  —L  +  jh,  the  quadrature  points  and  weights  are 


Hiql  =  hq- 1  +  h 
Hiq2  —  hq— i  -t-  h 

Hiq  3  =  hq- 1  +  h 


1  _  \J  15+2\/30 

2  2^35 


1  _  y/l5-2v/30 

2  2^35 


1  ,  \/l5-2v/30 

2 


2^35 


w  ,  =  _ Ml _ 

ql  12(18+ \/30) 


J  Wq2  = 


49/) 


12(18-t/30) 


(77) 


w  ,  = _ Mi _ 

’  qS  12(18-v/30) 


Hlq  4 


hq—  i  +  h 


1  .  \J  15+2\/30 

2  +  2^35 


Wqi  12(18+v/30) 


For  Nq  =  2  intervals,  and  hence  Nj  =  8,  the  quadrature  points  specified  by  (77)  are  depicted  in 
Figure  8(b). 

The  use  of  a  similar  relation  to  approximate  the  coercive  field  integral  yields  the  double  sum  (76) 
which  must  be  evaluated  when  computing  a  magnetization  value  for  each  input  field  value.  From 
(32),  it  is  observed  that  for  each  field  value  Hj.,  it  is  necessary  to  determine  whether  a  transition 
has  occurred  relative  to  the  coercive  value  HCi .  This  yields  IV*  x  Nj  conditions  to  be  checked  for  each 
input  value.  While  this  can  be  easily  accomplished  using  an  if-then  construct,  implementation  in 
this  manner  diminishes  significantly  the  efficiency  of  the  algorithm.  This  motivates  consideration  of 
an  algebraic  technique  for  evaluating  the  conditional  statements. 


4.5.2  Implementation  Algorithm  —  Hysteresis  Model  with  Negligible  After-Effects 

To  retain  the  history  of  whether  or  not  effective  field  values  Hjj  =  H  +  Hij  have  switched  due  to 
encounters  with  coercive  field  values  HCi ,  we  employ  (31)  to  motivate  the  matrix  formulation 

M=-  +  MrA(H;Hc,Hi)  (78) 

V 
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Figure  8:  (a)  Decay  exhibited  by  the  effective  field  Hj  having  the  density  e  Hi/2b 2  and  truncated 
domain  \—L,L\.  (b)  Gaussian  quadrature  points  •  and  initial  local  magnetization  values  for 

Nj  =  8.  (c)  Lognormal  density  u\ {Hc)  =  cie~^n(-Hc^Hc^2<^2  given  by  (64).  (d)  Distribution  of 

hysteresis  kernels  having  coercive  fields  Hc. 


where  A  =  1  if  evaluating  on  the  upper  branch  of  the  hysteron  and  A  =  —1  if  on  the  lower  branch. 
For  the  evaluation  of  (77),  A  is  an  TV,;  x  Nj  matrix  whose  ijth  component  specifies  whether  H j  has 
reached  the  coercive  value  HC/ .  We  also  define  the  following  matrices 


'  -l  •• 

•  -l 

i  •• 

•  i  ' 

tq 

1 _ 

_ 1 

^ init  — 

-l  •• 

•  -l 

i  •• 

■  •  i 

5 

Ni  x  Nj 

Hc  = 

_hCNi  ■ 

- 1 

sT 

... 

nk 


Hk  +  •  •  •  Hk  +  Hin,  I 


Hk  +  Hh 


Hk  +  HlN, 


Ni  x  Nj 


(79) 


and  weight  vectors 


WT  = 


-H2t  /2b2 
wie  /i/  ,  •••  ,wN.e 


-H2j2b2 


J  lx  AT,- 


ft  = 


Vie-[ln(Hcl/Hc)/2c]2 ,  ...  ’VNe-MHcN./Hc)/2c]* 


J  lxJVj 


(80) 


Here  Hk  =  H(tk )  is  the  kth  value  of  the  input  field.  The  magnetization  Mk  ~  M(Hk )  is  specified 
by  Algorithm  1.  Here  sgn  denotes  the  signurn  function,  .*  indicates  componentwise  multiplication, 

and  _  _ 

n  _  _  Cl  C2 

C  f  1  ' C2  h  bV2n 

follows  from  (71).  The  first  step  in  the  for- loop  updates  A  by  incorporating  the  status  of  previous 
coercive  field  switches.  We  note  that  by  employing  algebraic  matrix  operations  to  evaluate  and 
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incorporate  A,  the  efficiency  of  the  algorithm  is  improved  by  a  factor  of  more  that  100  over  algorithms 
utilizing  conditional  evaluations.  This  efficiency  is  crucial  for  system  design  and  real-time  control 
implement  ation . 

Algorithm  1. 

A  prev  —  A  init 

for  k  =  1  :  Nk 

A  —  Sgn(7"ffc  He-  *  A  prev) 

M  =  ^Hk  +  MrA 
Mk  =  CVfMW 

A prev  —  A 

end 

4.5.3  Implementation  Algorithm  —  Anhysteretic  Model 

The  implementation  of  the  discretized  anhysteretic  model  is  significantly  easier  than  implementation 
of  the  hysteresis  model  since  it  does  not  require  updating  of  A  to  retain  a  history  of  moment  switches 
due  to  local  coercive  fields.  One  can  employ  either  the  matrices  and  vectors  defined  in  (79)  and 
(80)  for  the  hysteresis  model,  or  a  reduced  set  of  vectors  which  reflects  the  fact  that  the  coercive 
density  integrates  to  the  constant  I\  defined  in  (71).  The  two  equivalent  approaches  are  illustrated 
in  Algorithms  2  and  3  where  C  =  I\C  and 

h-k  =  Hk  +  Hj  ,  •  •  •  ,  Hk  +  Hj  (82) 

1  vJl  xNj 

in  Algorithm  3.  Algorithm  2  retains  the  direct  correlation  with  the  hysteresis  model  whereas  Algo¬ 
rithm  3  is  more  efficient  to  implement  since  it  requires  vector  rather  than  matrix  multiplication. 

Algorithm  2. 

for  k  =  1  :  Nk 
A  =  sgn  (Hk) 

Man  =  ^k  +  MrA 
Mk  =  CVTManW 

end 

Algorithm  3. 

for  k  =  1  :  Nk 
A  =  sgn (hk) 

Man  =  ^hk  +  MrA 
Mank  =  CManW 

end 

5  Model  Properties  and  Validation 

To  illustrate  attributes  of  the  model,  we  consider  the  characterization  of  steel  under  zero  prestress 
conditions  using  data  from  [10].  This  demonstrates  a  variety  of  regimes  including  quantification 
of  the  anhysteretic  magnetization  and  the  prediction  of  symmetric  minor  loop  behavior  using  the 
model  with  parameters  obtained  from  the  symmetric  major  loop.  To  further  illustrate  the  manner 
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through  which  the  framework  characterizes  anhysteretic  material  behavior,  we  numerically  simulate 
the  experimental  process  in  which  a  decaying  AC  field  is  applied  at  a  DC  offset  to  demonstrate  that 
the  resulting  magnetization  coincides  with  that  predicted  by  the  anhysteretic  model.  In  the  third 
example,  we  illustrate  the  response  of  the  model  to  a  step  input  to  demonstrate  its  capability  for 
quantifying  the  effects  of  thermal  activation.  Finally,  we  note  that  additional  examples  illustrating 
the  capability  of  the  model  to  accurately  characterize  biased  minor  loop  behavior  can  be  found 
in  [63,64]  and  the  characterization  of  strains  in  a  Terfenol-D  transducer  is  illustrated  in  [36]. 

5.1  Experimental  Validation  for  Steel:  Anhysteretic  and  Biased  Minor  Loop 
Behavior 

We  consider  data  from  a  steel  specimen  having  a  length  of  6  cm  and  cross-sectional  area  of  1  cm  as 
reported  by  Jiles  and  Atherton  [10].  The  composition  (%  by  weight)  of  the  sample  was  C  (0.08), 
Mn  (1.98),  S  (0.08),  P  (0.015),  Cu  (0.055)  and  Mo  (0.235).  We  consider  hysteretic  and  anhysteretic 
data  collected  under  a  =  0  prestress  conditions  which  is  plotted  in  Figures  9(a)  and  9(b). 

Anhysteretic  Model 

The  anhysteretic  model  (72)  is  more  fundamental  than  the  full  hysteresis  model  (67),  since  it 
does  not  incorporate  local  coercive  fields,  so  we  consider  it  first.  The  parameters  Mr,  C.  q  and  b 
have  the  following  physical  interpretations.  The  local  remanence  value  Mr  and  constant  C  both 
scale  the  height  of  the  curve  and  are  constructed  to  yield  correct  saturation  values.  The  parameter 
r i  asymptotically  represents  the  reciprocal  slope  ^  at  field  reversal  and  an  initial  value  can  be 
obtained  from  the  slope  of  the  data  at  Hmax-  The  variance  b  quantifies  the  degree  of  pre-remanence 
switching  in  the  hysteresis  model  which  corresponds  to  the  curvature  of  the  anhysteretic  model  in 
low  drive  regimes. 

A  least  squares  fit  to  the  data  yielded  the  values  for  Mr,C,tj  and  b  summarized  in  Table  1.  A 
comparison  between  the  resulting  model  fit  and  data  in  Figure  9(a)  illustrates  that  the  anhysteretic 
model  (72)  accurately  characterizes  the  material  behavior  through  the  drive  range. 


(a)  (b) 


Figure  9:  Steel  data  from  [10]  collected  under  a  =  0  prestress  conditions,  (a)  Anhysteretic  data  and 
model  fit  provided  by  (72).  (b)  Hysteresis  data,  major  loop  model  fit  provided  by  (67),  and  minor 
loop  predictions. 
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Parameters 

V  Mr  (A/m) 

b  (A/m) 

C 

Hc  (A/m) 

c  (A/m) 

Values 

6.5  5.4  x  103 

3521.4 

0.0190 

250 

0.75 

Table  1:  Parameter  values  identified  for  the  anhysteretic  model  (72)  and  hysteresis  model  (67). 


Hysteresis  Model 

The  full  hysteresis  model  additionally  requires  the  estimation  of  the  mean  coercive  field  Hc  and 
variance  c.  Whereas  an  initial  estimate  for  the  former  can  be  obtained  directly  from  the  coercivity 
of  the  data,  the  parameter  c  which  quantifies  the  variability  at  coercivity  due  to  material  nonhomo¬ 
geneities,  is  usually  prescribed  a  qualitative  rather  than  quantitative  interpretation.  For  example, 
materials  exhibiting  a  steep  transition  at  coercivity  will  have  small  variances  when  compared  with 
materials  such  as  the  steel  data  plotted  in  Figure  9(b)  which  exhibits  a  gradual  transition  at  H  =  Hc. 

To  complete  the  model,  the  measured  coercive  field  Hc  =  910  A/m  was  used  as  an  initial  value  and 
the  values  of  Hc  and  c  compiled  in  Table  1  were  estimated  through  a  least  squares  fit  to  the  symmetric 
major  loop  data.  Measured  periodic  fields  having  lower  amplitudes  were  subsequently  input  to  the 
model  —  using  the  same  parameter  values  —  to  obtain  the  symmetric  minor  loop  predictions  which 
are  also  plotted  in  Figure  9(b).  It  is  observed  that  the  model  accurately  characterizes  the  hysteretic 
material  behavior  throughout  the  drive  regime,  including  the  approximately  quadratic  Rayleigh  loop 
behavior  at  low  input  fields. 


5.2  Numerical  Simulation  of  Anhysteretic  Behavior 

To  further  illustrate  the  manner  through  which  the  anhysteretic  magnetization  is  quantified  by 
this  framework,  we  numerically  simulate  the  experimental  procedure  used  to  obtain  Man  using  the 
full  hysteresis  model  (67),  and  compare  with  the  value  predicted  by  the  anhysteretic  model  (72). 
Specifically,  we  applied  the  periodic  and  subsequently  decaying  AC  field  depicted  in  Figure  10(a)  to 
the  model  (67)  to  simulate  the  experimental  procedure  used  to  obtain  Man  or  Ban  at  the  DC  field 
H0  =  2000  A/m.  The  parameter  values  from  Table  1  were  employed  so  the  result,  Ban  =  0.5544 
Tesla,  which  is  plotted  as  *  in  Figure  10(b),  is  representative  of  steel.  A  comparison  with  the 
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(a)  (b) 

Figure  10:  (a)  Field  input  to  the  hysteresis  model  (67)  to  generate  the  hysteretic  response  and 
anhysteretic  value  at  Ho  =  2000  A/m.  (b)  Value  of  Ban  generated  by  the  decaying  AC  field  (*)  and 
anhysteretic  model  (o),  and  full  anhysteretic  curve  ( - )  given  by  (72). 
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corresponding  value  of  Ban  =  0.5704  Tesla  predicted  by  the  anhysteretic  model  (72),  which  is  denoted 
by  o,  illustrates  that  the  two  approaches  yield  identical  results  to  within  reasonable  precision.  The 
locus  of  points  computed  using  (72)  completes  the  comparison  between  the  predicted  anhysteretic 
and  hysteretic  responses  for  the  material.  Hence  the  modeling  framework  developed  to  quantify 
hysteresis  in  ferromagnetic  materials  also  quantifies  the  anhysteretic  response  in  a  natural  manner. 

5.3  Experimental  Validation  for  Nickel:  Biased  Minor  and  Reversible  Post- 
Switching  Behavior 

To  illustrate  the  capability  of  the  model  to  characterize  biased  rnior  loop  behavior,  we  consider  data 
collected  from  a  rod  comprised  of  Nickel  200.  As  detailed  in  [65] ,  the  rod  had  a  diameter  of  0.0635  cm 
(1/4  in),  length  of  6.858  cm  (2.7  in),  and  was  employed  in  a  water-cooled  transducer  analogous  to 
the  Terfenol-D  design  depicted  in  Figure  2.  Data  was  collected  at  a  rate  of  10  mHz  to  minimize  eddy 
current  losses. 

The  kernel  (49),  which  incorporates  thermal  relaxation,  was  employed  in  the  magnetization  model 
and  general  densities  iq  and  zz2  were  identified  using  the  least  squares  techniques  detailed  in  [61]  for 
the  analogous  ferroelectric  model.  When  implementing  the  model  in  optimization  routines,  the 
measured  experimental  field  H,  plotted  in  Figure  11(a),  was  employed  as  input  to  the  model  to 
maintain  experimental  rates  and  time  scales.  The  resulting  magnetization  data  and  model  fit  are 
shown  in  Figure  ll(b)-(e).  The  time  history  in  Figure  11(b)  illustrates  the  fact  that  because  the 
nested  minor  loop  data  comprises  a  large  percentage  of  the  total  data  set,  the  optimization  routine 
determines  density  values  which  yield  greater  accuracy  in  the  minor  loops  than  in  the  post-switching 
region  of  the  major  loop.  The  slight  nonclosure  of  modeled  minor  loops  results  from  the  relaxation 
values  prescribed  during  the  optimization  process  to  accommodate  similar  behavior  in  the  measured 
data. 

For  fields  higher  in  magnitude  than  approximately  6  kA/rn,  the  model  incorporates  the  reversible 
post-switching  behavior  of  the  material  due  to  the  form  of  the  energy-based  kernel.  The  characteri¬ 
zation  of  these  effects  using  a  Preisach  approach  requires  extensions  of  the  type  detailed  in  [9]  since 
classical  Preisach  hysterons  have  zero  slope.  Similarly,  an  extended  formulations  of  the  Jiles- Atherton 
model  is  required  to  incorporate  this  effect. 

The  accuracy  exhibited  in  Figure  11(b)  is  important  for  model-based  control  design  since  charac¬ 
terization  and  compensation  in  this  context  are  typically  posed  as  a  functions  of  time.  Further  details 
regarding  the  utility  of  the  framework  for  model-based  control  design  are  provided  in  Section  6. 

5.4  Numerical  Simulation  of  Relaxation  Behavior 

One  manifestation  of  relaxation  effects  is  the  nonclosure  of  biased  minor  loops  observed  in  Figure  12. 
To  further  illustrate  the  manner  through  which  the  model  characterizes  magnetic  after-effects,  the 
discontinuous  step  input  shown  in  Figure  12(a)  was  input  to  the  model.  The  resulting  magnetization, 
plotted  in  Figure  12(b),  exhibits  both  an  associated  discontinuity  and  creep  or  after-effects  due  to 
the  modeled  thermal  activation  mechanisms.  Quantification  of  this  phenomenon  is  important  both 
for  fundamental  material  modeling  and  characterization  of  transducers  required  to  hold  a  precise  set 
point  over  time  scales  comparable  with  relaxation  times. 

For  the  classes  of  macroscopic  models  discussed  in  Section  1,  extended  Preisach  formulations  in¬ 
corporating  Arrhenius  behavior  will  characterize  this  behavior  whereas  the  Jiles-Atherton,  Globus, 
and  Stoner-Wohlfarth  formulations  do  not  presently  address  this  phenomenon.  The  inherent  in¬ 
corporation  of  thermal  activation  mechanisms  constitutes  an  advantage  of  the  homogenized  energy 
framework. 
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Figure  11:  (a)  Field  input,  (b)  Nickel  200 
H-M  data,  (d)  H-M  model  fit,  and  (e)  con: 
following  negative  remanence. 
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Figure  12:  (a)  Input  field  H  and  (b)  response  of  the  model  (69)  employing  the  kernel  (49)  which 
incorporates  thermal  activation. 


6  Concluding  Remarks 

The  model  presented  here  provides  a  framework  for  characterizing  hysteretic  and  anhysteretic  be¬ 
havior  of  a  broad  range  of  ferromagnetic  compounds.  The  approach  combines  energy  analysis  at 
the  lattice  level  with  stochastic  homogenization  techniques  to  construct  macroscopic  models  which 
accurately  characterize  the  nonlinear  H-M  or  H-B  material  behavior  while  remaining  sufficiently 
efficient  to  permit  subsequent  transducer  design  or  model-based  control  designs  having  the  potential 
for  real-time  implementation.  At  the  lattice  level,  exchange  and  magnetostatic  energy  relations  are 
employed  in  a  mean  field  framework  to  construct  a  temperature-dependent  Gibbs  energy  relation 
G  which  quantifies  the  transition  from  ferromagnetic  to  paramagnetic  behavior  as  temperatures  are 
increased  through  the  Curie  temperature  Tc.  For  fixed  temperature  regimes,  Taylor  expansion  about 
equilibria  yields  a  piecewise  quadratic  Gibbs  relation  which  is  highly  efficient  to  implement.  For 
regimes  in  which  magnetic  after-effects  (thermal  relaxation)  are  negligible,  the  local  average  magne¬ 
tization  M  at  the  lattice  level  is  determined  from  the  necessary  condition  =  0  whereas  Boltzmann 
principles  are  employed  to  balance  G  with  the  relative  thermal  energy  kT /V  for  operating  regimes 
in  which  relaxation  times  are  significant  compared  with  drive  frequencies. 

To  construct  macroscopic  models  utilizing  the  energy-based  kernels  or  hysterons  M,  we  assume 
that  bias  fields  Hh ,  local  coercive  fields  HCl  and  interaction  fields  Hj  are  manifestations  of  underly¬ 
ing  distributions  rather  than  constant  parameters  or  inputs.  This  yields  constitutive  relations  which 
guarantee  minor  loop  closure  once  accommodation  and  after-effects  are  completed,  incorporate  relax¬ 
ation  and  certain  temperature-dependencies,  incorporate  reversible  effects,  and  yield  noncongruency 
in  certain  operating  regimes.  The  model  is  also  highly  efficient  to  implement  thus  facilitating  future 
incorporation  in  design  and  control  algorithms.  Finally,  the  framework  incorporates  mechanisms 
common  to  a  number  of  presently  employed  macroscopic  models. 

J lies- Atherton:  The  Ising  relation  (29)  obtained  by  minimizing  the  Gibbs  energy  constructed 
through  statistical  mechanics  tenets  is  analogous,  and  agrees  through  first-order  terms,  with  the 
Langevin  relation  (1)  employed  in  the  Jiles- Atherton  theory  to  quantify  the  anhysteretic  magneti¬ 
zation.  The  difference  between  the  theories  lies  in  the  manner  through  which  energy  relations  are 
constructed  and  losses  are  incorporated. 

Globus  and  Stoner-Wohlfarth:  The  energy-based  hysterons  M  provided  by  the  proposed  model 
are  also  analogous  to  those  provided  by  the  Globus  and  Stoner-Wohlfarth  models.  The  present  theory 
differs  from  the  original  Globus  and  Stoner-Wohlfarth  theories  due  to  its  formulation  in  terms  of  the 
exchange  energy  which  incorporates  moment  interactions.  It  also  differs  from  the  Stoner-Wohlfarth 
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theory  in  that  it  does  not  incorporate  anisotropy  energies. 

Preisach:  The  theory  bears  the  greatest  resemblance  to  extended  Preisach  theories  and  the 
formulation  (69)  in  terms  of  general  densities  provides  an  energy  basis  for  certain  extended  Preisach 
models.  The  primary  differences  between  the  proposed  framework  and  classical  Preisach  theory  are 
the  following,  (i)  As  detailed  in  Section  5.1,  certain  parameters  such  as  the  local  average  coercive 
field  Hc  and  reciprocal  slope  i]  at  Hmax  can  be  directly  correlated  with  properties  of  the  data  to 
facilitate  model  construction  and  updating.  This  is  analogous  to  interpretations  proposed  in  [9]  for 
Preisach  models  with  a  priori  density  specifications,  (ii)  Temperature  and  relaxation  mechanisms 
are  incorporated  in  the  basis  or  kernel  M  of  (67)  or  (69)  rather  than  the  densities  or  weights  as  is 
typically  the  case  for  Preisach  formulations.  This  eliminates  the  need  for  vector- valued  measures 
or  lookup  tables  for  general  operating  conditions,  (iii)  As  illustrated  in  Figure  5,  the  proposed 
framework  employing  the  kernel  (29)  relaxes  the  criterion  of  congruency  and  hence  is  analogous  to 
certain  moving  Preisach  models  [9] .  (iv)  The  model  also  automatically  incorporates  reversible  effects 
after  switching  and  hence  does  not  required  the  Preisach  extensions  detailed  in  [9]  to  incorporate 
reversibility. 

We  note  that  the  framework  developed  here  for  ferromagnetic  materials  is  based  on  theory  quan¬ 
tifying  phase  transitions  in  shape  memory  alloys  (SMA)  [66-69]  which  has  subsequently  been  ex¬ 
tended  to  ferroelectric  materials  [61,70].  Hence  the  framework  provides  a  unified  methodology  for 
modeling  constitutive  nonlinearities  and  hysteresis  in  ferroelectric,  ferromagnetic  and  ferroelastic 
materials  [63,64].  This  directly  facilitates  the  construction  of  unified  characterization  techniques 
for  ferroic  compounds  and  the  development  of  unified  model-based  control  designs  for  ferroic  trans¬ 
ducers.  It  also  provides  a  framework  which  may  facilitate  the  construction  of  models  for  hybrid 
transducers  employing  multiple  components  (e.g.,  Terfenol-D  and  PZT)  or  the  characterization  of 
newly  developed  and  proposed  materials  such  as  ferromagnetic  shape  memory  alloys  (FSMA)  which 
exhibit  attributes  common  to  multiple  classes  of  the  compounds. 

As  noted  in  Section  2,  exchange,  magnetostatic,  magnetoelastic  and  anisotropy  energy  compo¬ 
nents  contribute  to  the  total  energy  in  ferromagnetic  materials.  The  present  modeling  framework 
incorporates  the  exchange  and  magnetostatic  energies,  and  aspects  of  the  magnetoelastic  energy  have 
been  employed  in  [36]  to  develop  magnetoelastic  constitutive  relations 

a  =  YMe  -  Ym^M2 

POO  POO  (83) 

M(H)=  /  /  viHcH^MiH  +  H^Hc^dHjdHc 

J  0  J  —  oo 

where  e,  a  and  7  respectively  denote  a  uniaxial  strain,  stress  and  coupling  coefficient.  These  relations 
form  the  basis  for  constructing  PDE  models  for  magnetostrictive  transducers,  and  the  performance  of 
the  framework  for  characterizing  the  magnetization  and  strains  generated  by  a  Terfenol-D  transducer 
analogous  to  that  depicted  in  Figure  2  is  demonstrated  in  [36]. 

Attributes  and  open  research  questions  pertaining  to  the  framework  can  be  summarized  as  fol¬ 
lows.  (i)  Minor  Loop  Closure  and  After-Effects:  For  operating  regimes  in  which  thermal  relaxation 
processes  are  negligible,  use  of  the  kernel  (32)  in  the  macroscopic  model  (69)  guarantees  the  closure 
of  biased  minor  loops.  Alternatively,  use  of  the  kernel  (49),  which  incorporates  thermal  relaxation, 
yields  nonclosure  of  biased  minor  loops  in  certain  regimes  as  well  as  after-effects  of  the  type  shown  in 
Figure  12.  (ii)  Noncongruency:  As  noted  in  Remark  2  and  depicted  in  Figure  5(a),  the  model  employ¬ 
ing  the  kernel  (29)  incorporates  the  noncongruency  exhibited  by  certain  materials,  (iii)  Reversibility: 
The  model  employing  the  kernels  (32)  or  (49)  provides  reversible  behavior  after  switching  due  to  the 
nonzero  slope  of  the  upper  and  lower  hysteron  branches  (e.g.,  see  Figure  11).  (iv)  Accommodation: 
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Whereas  the  framework  phenomenologically  quantifies  certain  accommodation  processes,  the  exten¬ 
sion  of  the  theory  to  incorporate  energy  mechanisms  associated  with  magnetic  accommodation  has 
not  been  completed  and  is  under  present  investigation,  (v)  Temperature  and  Stress- Dependencies: 
The  kernel  (29)  incorporates  certain  temperature-dependencies  as  well  as  the  transition  between 
ferromagnetic  and  paramagnetic  phases;  however,  its  accuracy  should  be  considered  limited  when 
quantifying  changes  over  a  broad  temperature  range.  Certain  extensions  to  the  framework  to  provide 
more  comprehensive  characterization  of  temperature  effects  has  been  developed  in  [35]  in  the  context 
of  relaxor  ferroelectric  materials.  The  framework  presented  here  focuses  on  constant  stress  condi¬ 
tions.  Initial  extensions  of  the  theory  to  incorporate  stress-dependencies  in  M  due  to  magnetoelastic 
coupling  are  provide  in  [73].  (vi)  Eddy  Currents:  The  present  framework  does  not  incorporate  eddy 
currents  and  hence  should  be  employed  in  low  frequency  regimes  or  transducers  constructed  for  mini¬ 
mal  eddy  losses  (e.g.,  laminates),  (vii)  Anisotropy:  The  model  is  presently  formulated  for  isotropic  or 
weakly  anisotropic  materials  and  the  extension  of  the  framework  to  incorporate  the  energy  associated 
with  hexagonal  and  cubic  crystalline  anistropies  is  under  investigation. 

One  of  the  motivating  criteria  when  constructing  the  characterization  framework  was  to  provide 
macroscopic  models  with  sufficient  efficiency  to  permit  model-based  control  design.  It  is  illustrated  in 
[71]  that  the  nronotonicity  inherent  to  the  H-M  relation  can  be  exploited  to  construct  highly  efficient 
approximate  inverse  representations  M_1(i7)  using  the  discretized  model  (76).  These  approximate 
inverses  are  then  employed  as  filters  to  the  linear  and  hysteretic  transducers  so  that  prescribed 
and  actual  control  inputs  to  the  plant  approximately  coincide.  It  is  illustrated  in  [71,  72]  that 
by  linearizing  a  Terfenol-D  transducer  in  this  manner,  highly  accurate  tracking  tolerances  can  be 
maintained  through  linear  robust  control  designs  while  operating  in  highly  nonlinear  and  hysteretic 
drive  regimes. 

7  Appendix  A 

We  summarize  here  the  attributes  of  a  Dirac  sequence  and  provide  a  theorem  which  establishes  that 
the  convolution  of  Dirac  sequences  with  suitably  smooth  functions  /  will  limit  to  a  point  evaluation 
of  /.  This  theorem  is  employed  in  the  model  development  to  illustrate  the  convergence  of  Gaussian 
moment  distributions  to  a  Dirac  distribution  in  the  limit  of  small  relative  thermal  energies  kT /V  or 
increasing  control  volumes  V.  It  also  provides  a  framework  employed  to  demonstrate  that  models 
which  include  thermal  after-effects  converge  to  models  based  on  the  assumption  of  negligible  thermal 
energy  as  control  volumes  are  taken  to  be  arbitrarily  large. 

Theorem  1.  Let  {4>j}  be  a  sequence  satisfying  the  following  properties: 

(i)  4>j  >  0  for  all  j 

(ii)  f  <j>j(y)dy  =  1  for  all  j 

(Hi)  Given  e,5  >  0,  there  exists  jo  such  that 


fij(y)dy  <  e 


for  all  j  >  jo. 

Let  f  be  piecewise  continuous  on  M,  continuous  on  the  interval  [a,  b] ,  and  satisfy  the  decay  property: 
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(iv)  Given  s,  S  >  0,  there  exists  jo  such  that 


for  all  j  >  jo  and  x  6  [a,  b] . 


My)f(x 


y)dy  <  e 


Then  for  x  6  [a,  b\,  f>k  *  f  converges  to  f;  that  is 


J  <Pj(y)f{x  -  y)dy  -»•  f(x) 

or 

J  </>j(x-y)f(y)dy^f(x). 


Proof.  From  (ii)  it  follows  that 

f(x)  =  f(x)  J  <f>j(y)dy  =  J  f(x)<t>j(y)dy 

so  that  for  x  €  [a,  b] , 

<t>j*  f(x)  -  f(x)  =  J  f>j{y)f{x-y)dy-  J  </>j(y)f(x)dy 
=  J  &j{y)[f(x-y)  -f(x)]dy. 

From  the  continuity  of  /,  it  follows  that  for  fixed  e,  there  exists  5  such  that 

I  f(x-y)  ~  f(x)  |  <  £ 

for  \y\  <5.  For  this  5,  we  write 


\4>j  *  f(x)  -  f(x)\  <  [  </>j(y)\f(x  -y)  -  f(x)\dy  + 

/  4>j{y)f{x  -  y)dy 

J\y\<5 

J\y\>s 

+  /  f>j{y)\f{x)\dy. 

J\y\>8 

For  sufficiently  large  j,  the  integral  over  the  region  |y|  <  5  is  bounded  by  e  due  to  the  continuity 
of  /  on  [a,  b]  whereas  the  second  integral  is  bounded  by  e  due  to  (iv).  Finally  the  third  integral 
is  bounded  by  ||/||oo£,  where  ||/||oo  =  maxxeraM  \f(x)\,  thus  yielding  the  desired  convergence.  The 
convolution  expression  follows  from  a  direct  change  of  variables.  □ 


We  note  that  a  sequence  of  functions  {4>j}  satisfying  the  properties  (i)-(iii)  is  termed  a  Dirac  sequence 
on  M1.  Additionally,  if  we  replace  the  assumption  (iv)  by  the  condition  that  /  is  bounded  and 
measurable  on  R,  then  Theorem  1  is  a  1-D  version  of  Theorem  3.1  from  page  228  of  [74]. 
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